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SUMMARY 


The primary objective of this research program is to improve our 
understanding of the physical processes occurring in the interaction of the solar 
wind with Venus. This will be accomplished through the use of numerical solutions 
of the two- and three-dimensional magnetohydrodynamic (MHD) equations and 

through comparisons of the computed results with available observations. A large 
portion of this effort involves the study of processes due to the presence of the 
magnetic field and the effects of mass loading. 

Data from the Pioneer Venus Orbiter (PVO), as well as data from Venera 9 and 
10 and other sources, have contributed to a growing understanding of the interaction 

of the solar wind with an unmagnetized body such as Venus. The general picture 
derived from the available data is that, for typical solar wind conditions, an 

electrically conducting ionosphere deflects the supersonic solar wind around the 

planet. A bow shock forms upstream of the ionosphere and serves to slow, and also 

assists in deflecting, the solar wind. The boundary separating the chocked and 
magnetized solar wind plasma from the thermal ionosphere plasma is referred to as 
the ionopause, and the region between the bow shock and the ionopause is referred 
as the magnetosheath. The inner portion of the magnetosheath contains a region of 
enhanced magnetic pressure referred to as the magnetic barrier. It is the physical 

processes involved in these relatively large-scale interaction region as well as in the 
downstream magnetic tail that this work has application to. 

The general approach that we use is to treat the planet as a solid, ideally 
conducting electrical sphere. Consequently, the plasma velocity and magnetic field 

are required to be parallel to the surface of the planet. The equations are solved in 

time and space in a spherical coordinate system centered on the planet. A particular 

simulation is obtained by maintaining predetermined, fixed solar wind conditions on 
the dayside of the outer boundary. For given initial conditions the computation is 
continued in time until the interaction of the specified solar wind with the planet 

has reached an equilibrium in time. The spatial distribution of the initial velocity 
and magnetic field are taken from the potential solutions for incompressible flow 
over a sphere. 

1. Scope of the investigation 

The four specific research topics that have been investigated in detail to date 
include (1) the study of bow shocks at unusually large distances from the planet, (2) 

the effect of the interplanetary magnetic field direction and strength on the 

interaction, (3) the effect of mass loading on the interaction for the special 
orientation of the interplanetary magnetic field perpendicular to the flow velocity, 
and (4) the effect of mass loading for other orientations relative to the flow velocity 
and for other values of the plasma beta. Each of the four studies carried out are 
described in more detail below. All of these results have been presented in 

international or national meetings and published in journals. 

2. Venus bow shocks at unusually large distances from Venus 

Recent analysis of data from the Pioneer Venus Orbiter (PVO) has shown that 


1 




the bow shock often travels to unusually large distances from the planet when the 
solar wind magnetosonic Mach number is near unity. We suggest that distant bow 

shocks can be explained as an integral part of the response of the global solar 

wind/Venus interaction to anomalous local solar wind conditions that existed during 

the time of these observations. The lower than normal plasma beta and 
magnetosonic Mach number are in a parameter regime for which the usual fast- 

mode bow shock close to the planet may not provide the necessary compression and 
deflection of the solar wind. The solar wind conditions in the vicinity of Venus 
become such that the usual fast MHD shock near the planet is no longer a stable 

solution. The fast MHD shock is replaced by a shock configuration containing an 

intermediate MHD shock over at least a portion of the shock front, and the resulting 

configuration propagates upstream. 

The MHD equations without explicit dissipation are used in this study. The 

equations are solved numerically using the modified Lax-Wendroff scheme with a 

smoothing term. In our simulations the planet is replaced by a cylinder aligned 

along the z-axis in a cylindrical coordinate system. The z-axis is taken to be in the 
ecliptic plane perpendicular to the Sun-Venus plane, and the MHD equations are 
solved in the r-9 plane. The cylindrical approximation should retain the essential 

physics for this problem. 

Using MHD simulations we show that, for these conditions, the usual fast shock 
is replaced by a bow shock consisting of an intermediate shock near the Sun-Venus 

line and a fast shock at large distances from the Sun-Venus line. This composite bow 

shock propagates upstream away from the planet at a low speed and appears to be 

approaching a new equilibrium stand-off location at a large distance from the 

planet. These results are published in the article entitled “Venus bow shocks at 

unusually large distances from the planet” in Geophysical Research Letter, Vol 20 

755-758, 1993. 

3. Three-dimensional MHD simulations of the interaction between Venus and the 

solar wind 

We have developed a three-dimensional, single-fluid, MHD code which has 
successfully simulated the major qualitative features of the Venus-solar wind 
interaction. We model the Venus obstacle as a hard, conducting sphere and the solar 
wind as a single-fluid MHD plasma. We simulate the case in which the IMF is tilted 
with respect to the incoming solar wind velocity at approximately the Parker spiral 
angle; this is the first simulation of this particular case which includes all spatial 
regions of the interaction. A realistic value of beta has also been used. 

Our simulations show that a bow shock forms at a height above the ionosphere 
close to the that deduced from Pioneer Venus Orbiter (PVO) observations. A magnetic 
barrier beneath the bow shock is found. The bow shock shows an asymmetry in the 
terminator plane qualitatively similar to results obtained from PVO data analyses. 
The magnetic field drapes over the planet and forms a two-lobed magnetotail. The 
magnetic field configuration near the planet is very similar to that recorded by PVO 
as well. 


Our simulations have also produced several interesting results which seem to 
have not been studied before. In the asymmetric 45° simulation the peak in the 
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magnetic field buildup is shifted from the subsolar point toward the quasi- 
perpendicular region of the shock. In contrast, the peaks in thermal pressure and 
density are shifted toward the quasi-parallel region. In the near-tail region our 
perpendicular simulation has established the possibility of magnetic reconnection. 


These results are published in the article entitled “Three-dimensional MHD 
simulations of the interaction between Venus and the solar wind” in Journal of 
Geophysical Research, vol. 100, 21645-21658, 1995. 

4. 2-D Numerical simulations of mass loading in the solar wind interaction with 

Venus 

Numerical simulations are performed in the frame work of nonlinear two- 
dimensional magnetohydrodynamics to investigate the influence of mass loading on 
the solar wind interaction with Venus. A modified Lax-Wendroff scheme is used to 
solve the equations in a spherical coordinate system. For these computations the 
detailed chemistry of the ionosphere is neglected, and the planet is treated as a 
conducting sphere. Numerical computations have been performed to study the effect 
of mass loading and the IMF strength on the magnetic barrier and the general 
configuration of the magnetic tail. The IMF orientation included in this study is 

parallel to the solar wind flow. 

The principal physical features of the interaction of the solar wind with the 

atmosphere of Venus are presented. The formation of the bow shock, the magnetic 
barrier, and the magnetotail are some typical features of the interaction. The 
deceleration of the solar wind due to the mass loading near Venus is an additional 
feature. The effect of the mass loading is to push the shock farther outward from the 
planet. The influence of different values of the magnetic field strength on the 
plasma evolution is considered. 

The main results are the following: The mass loading increases the mass 

density in the overall region, except the terminator and tail where the mass density 

depletion occurs as a consequence of the nonlinearity action. The mass density 
depletion is larger for a higher beta. As a consequence of this reduction, the plasma 
flow, pressure, and temperature are reduced. This scenario depends, however, on the 
value of the plasma beta. It was found that for beta = 0.1 both the plasma pressure 
and the temperature are increased at the dayside by the mass loading. 

These results were published in the article entitled “Numerical simulations of 
mass loading in the solar wind interaction with Venus” in Journal of Geophysical 
Research, Vol. 101, 2547-2560, 1996. 

5. 3-D Numerical modeling of the solar wind interaction with Venus 

The large scale (global) interaction of atmosphere of Venus with the 
surrounding medium is numerically simulated in the framework of 3 -dimensional 
magnetohydrodynamics. This study extends the analysis described in Section 4 into 
the case of the IMF perpendicular to the solar wind flow. The numerical simulations 
have been carried on for a long time which allow us to believe that the simulation 
approached a quasi-steady state. 
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The simulations are performed for the case of the interplanetary magnetic 
field perpendicular to the solar wind flow. The model reproduces several features of 

the interaction predicted by earlier theories and observations. The solar wind 
interaction with Venus is characterized by the mass loading of the solar wind with 
heavy oxygen ions which are produced by the photoionization of neutrals in the 
extensive ionosphere. This mass loading slows down the solar wind and ultimately 
leads to the formation of a magnetic barrier and a magnetotail. The location of the 
bow shock is barely affected by different values of the magnetic field strength. The 
numerical results indicate that the presence of mass loading causes the bow shock to 
move outward from the planet but not so far as observed at Venus during the period 
of maximum solar activity. 

The perpendicular case, although similar in some aspects, differs in a general 
behavior to the parallel case. In particular, it is found that the plasma dynamics is 

hardly affected by the strength of the magnetic field which originated from the 
solar wind. The location of the bow shock is virtually the same for different values 

of the plasma beta. The IMF lines pile up against the atmospheric plasma and drape 
over the planet to form a magnetotail with two lobes of opposite magnetic polarity 
separated by a plasma sheet. The parallel magnetic field causes the plasma dynamics 
to be highly sensitive to the value of the plasma beta. 

These results were published in the article entitled “Numerical modeling of 
the solar wind interaction with Venus” in Planetary Space Science, Vol. 44, 243-252 
1996. ‘ 
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VENUS BOW SHOCKS AT UNUSUALLY LARGE DISTANCES FROM THE PLANET 

R. S. Steinolfson and S. Cable 

Department of Space Sciences. Southwest Research Institute, San Antonio, Texas 


Abstract . Recent analysis of data from the Pioneer Venus 
Orbiter (PVO) has shown that the bow shock often travels to 
unusually large distances from the planet when the solar wind 
magnetosonic Mach number is near unity. We suggest that 
distant bow shocks can be explained as an integral part of the 
response of the global solar wind/Venus interaction to the 
anomalous local solar wind conditions that existed during the 
time of these observations. The lower than normal plasma 
beta and magnetosonic Mach number are in a parameter 
regime for which the usual fast-mode bow shock close to the 
planet may not provide the necessary compression and 
deflection of the solar wind. Using MHD simulations we 
show that, for these conditions, the usual fast shock is 
replaced by a bow shock consisting of an intermediate shock 
near the Sun-Venus line and a fast shock at large distances 
from the Sun-Venus line. This composite bow shock propa- 
gates upstream away from the planet at a low speed and 
appears to be approaching a new equilibrium stand-off 
location at a large distance from the planet 

Introduction 

For typical solar wind conditions the Venus bow shock is 
located within or near one radii (R v ) of the planet surface 
along the Sun-Venus line. However, for low magnetosonic 
Mach numbers during two orbits of PVO the bow shock was 
observed out to approximately 12 R v away fron the planet at 
a 45° angle from the Sun-Venus line for one orbit and at 11 
R y above the terminator for the other [Russell and Zhang, 
1992]. These low magnetosonic Mach numbers were due 
primarily to an enhanced interplanetary magnetic field (about 
a factor of 3 larger than more typical values) and also to a 
reduced solar wind density (about a factor of 3 lower than 
normal). Neither the solar wind speed nor the temperature 
deviated appreciably from average values; hence, the sound 
Mach number remained near a typical value. As a result, the 
plasma beta was lower than normal and most likely became 
less than or even much less than unity. The plasma beta can 
be written in terms of the Alfv£n Mach number M b and the 
sound Mach number M s as P=2(M b /M J ) 2 /y. The poly tropic 
index (ratio of specific heats) y is generally taken as 5/3. In 
a low-beta environment (p < l) the fast-mode magnetosonic 
Mach number becomes equal to the Alfvdn Mach number for 
parallel propagating shocks. 

Russell and Zhang [1992] estimated that the Alfvdn and 
magnetosonic Mach numbers were unity to within their 
ability to measure them. For a temperature of 10 5 °K and a 
solar wind speed of 500 km/sec (typical of conditions during 
their observations), the sound Mach number becomes approxi- 
mately 9.5. Consequently, the plasma beta is about 0.013. For 
more usual solar wind conditions near Venus beta would have 
a value near 1.5. 

We suggest that a possible explanation for the above 
erratic behavior of the bow shock is that the solar wind 
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conditions in the vicinity of Venus become such that the 
usual fast MHD shock near the planet is no longer a stable 
solution. The fast MHD shock is replaced by a shock 
configuration containing an intermediate MHD shock over at 
least a portion of the shock front, and the resulting configura- 
tion propagates upstream. 

Formation of Intermediate Shocks in CMEs 

We begin by noting the similarity between bow shocks 
formed in the interaction of the solar wind with Venus and 
shocks formed ahead of CMEs in the solar corona (illustrated 
for the formation of fast shocks in the schematic in Figure 1). 
The analogy is most apparent when the north-south compo- 
nent of the interplanetary magnetic field (IMF) is small as 
assumed in the figure. The bow shock under normal condi- 
tions is approximately stationary in the laboratory frame, 
while the shock ahead of a CME propagates out through an 
expanding plasma. The similar behavior for the two examples 
of shock formation is most apparent in the rest frame of the 
individual shocks. A spiral component for the IMF (perpen- 
dicular to the figure) would not dramatically influence the 
analogous behavior for the two situations. Clearly, the 
essential physics should be similar for the two cases even for 
small north-south or spiral components of the IMF. 

For given values of the shock speed and upstream condi- 
tions, the MHD shock-jump (Rankine-Hugoniot) equations 
reduce to a cubic polynomial and formally allow three 
physically realistic shock solutions with an entropy rise (slow, 
intermediate and fast shocks). It was long believed that 
intermediate shocks did not exist in nature and were nonevo- 
lutionary or extraneous [e.g., Kantrowitz and Petschek, 1966]. 



Fig. 1. A schematic illustrating the similarity between (a) the 
formation of a fast MHD bow shock in the interaction of the 
solar wind with Venus and (b) the formation of a fast MHD 
shock ahead of a CME propagating out through the solar 
corona. 
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Recent studies, however, iend support to the possibility that 
ihey may not be extraneous [e.g„ Wu.4990]. Steinolfson anu 
Hundhausen [1990] show that intermediate shocks are 
necessary to provide a smooth and continuous transition from 
fast shocks to slow shocks with decreasing CME speec. 

An important difference between the three types of MHD 
shocks is the etfect that the shock has on the component of 
the magnetic field perpendicular to the shock normal. This 
component maintains the same direction across both slow and 
last shocks. It increases in magnitude across a fast shock and 
decreases in magnitude across a slow shock. That is, fast 
shocks bend the magnetic field toward the shock plane, and 
slow shocks bend the magnetic field away from the shock 
plane. The perpendicular component of the magnetic field 
reverses direction across an intermediate shock, and the field 
magnitude may either increase or decrease. 

In order for a particular type of shock to form, the shock 
speed (or solar wind speed in the case of bow shocks) must 
exceed the respective linear wave speed upstream of the 
shock. The wave speeds depend on the angle 0 between the 
shock normal and the upstream magnetic field. For the 
relatively idealized configurations shown in Figure I, the 
relevant wave speeds are close to the values for parallel 
shock propagation. The slow speed reduces to the sound 
speed at parallel propagation, and the fast speed becomes the 
Altvdn speed. The intermediate speed also reduces to the 
Altvdn speed at parallel propagation when pel. Consequently, 
lor the physical conditions we are interested in, intermediate 
shocks only form when the Alfvgn Mach number is larger 
than 1 (M b >l). 

As shown by Steinolfson and Hundhausen [1990], there is 
also a maximum Alfv6n Mach number for which intermediate 
shocks may form given by M\ cnl = [7(l-p)+l]/(Y-l) t where 
P is the value upstream of the shock. This expression was 
derived by Kennel and Edmiston [1988] as the maximum 
value for the formation of switch-on shocks. However, the 
formation of switch-on shocks and the formation of interme- 
diate shocks are intimately related [Steinolfson and Hund- 
hausen, 1990], As noted in the introduction, for small values 
of P (=10 MO 1 ) the sound Mach number may be large even 
though the Alfvdn Mach number is relatively small (near 
unity). The maximum value of M bcm occurs at P=0, in which 
case it has a value of 2 for y=5/3, so the shock speed cannot 
be too high. Intermediate shocks do not exist at ail for P>1.2. 
The net result is that intermediate shocks form when the 
CME shock Alfv£n Mach number is in the range 1 < M b < 
Mbcnt- By analogy, the solar wind Alfvtfn Mach number must 
be in the above range for intermediate shocks to form along 
a portion ot the bow shock. This is exactly the parametric 
regime of the solar wind parameters for both P and M b when 
the distant bow shocks were observed at Venus. 

When intermediate shocks do form near the leading edge 
ot CMEs, the shock Iront has the configuration shown in the 
schematic in Figure 2 [Steinolfson and" Hundhausen. 1990]. 
The shock is assumed to be propagating vertically along the 
vertical magnetic field in the tigure. The central portion of 
the shock tront that is bowed downward away from the 
direction ot travel (between the two switch-on shocks) 
contains the intermediate shocks. The shocks at the Hanks 
beyond the switch-on shocks are fast shocks, and a gas 
dynamic shock exists at the point of parallel propagation at 
the center of the shock front. 

MHD Simulation of Distant Bow Shocks 

The MHD equations without explicit dissipation are used 
:n this study. The equations are solved numerically using the 



Fig. 2. A schematic illustrating the shape of a CME shock 
front when an intermediate MHD shock forms along a portion 
of the shock front. The lines with an arrow at the top indicate 
magnetic Field lines, and the parallel curved lines are the 
shock front. The shock is travelling vertically along the 
uniform and parallel magnetic field lines. Note that the 
component of the magnetic field perpendicular to the shock 
normal reverses direction across the intermediate shock 
portion of the shock front located between the two switch-on 
shocks at parallel propagation. 


modified Lax-Wendroff scheme given by Rubin and Burstein 
[1967] with a smoothing term suggested by Lapidus [1967]. 
In our simulations the planet is replaced by a cylinder aligned 
along the z-axis in a cylindrical coordinate system. The z-axis 
is taken to be in the ecliptic plane perpendicular to the Sun- 
Venus line, and the MHD equations are solved in the r,0 
plane. This cylindrical approximation should retain the 
essential physics for this problem, although more realistic 
simulations would, of course, use a three-dimensional 
geometry. 

The simulation is carried out for a cross-section of a 
cylindrical annulus lying between an inner radial boundary at 
the surface of the planet R v and an outer radial boundary 
taken to be at 15 R,. The 0 axis is selected so that 0=0° is 
directed toward the Sun along the Sun- Venus line. The solar 
wind velocity and fMF are chosen such that the physical 
quantities are either symmetric (for p, p, B r ) or anti- 
symmetric (for u<j, B e ) about the 0=0° and the 0=180° axes so 
the computation must be computed only for the numerical 
box given by Rv < r < 15 R* and 0° < 0 £ 180°. The angular 
grid spacing is 1.5°, and the radial spacing is 0.15 R^ giving 
a grid of 122 by 95. The plasma velocity and magnetic field 
are required to be parallel to the surface of the planet, and 
zero-order radial extrapolation is used to determine the 
thermodynamic quantities on the planet surface. The physics 
and chemistry of the ionosphere should not be significant for 
the present study and are not included. Solar wind conditions 
are maintained on the dayside of the outer radial boundary, 
and zero-order extrapolation along the local flow direction is 
used to update values on the nightside. 

The simulation is initiated without the magnetic field with 
the solar wind thermodynamic conditions specified throughout 
the numerical box and the solar wind velocity specified 
everywhere beyond r=2R v . When this gas dynamic computa- 
tion reaches a dynamic equilibrium, the initial magnetic field, 
as given by the potential flow solution for flow over a 
cylinder, is specified at every grid point. The time-dependent 
computation is then restarted and follows the solution as far 
in time as desired. 

The results of two separate simulations are shown below. 
For the first one the solar wind conditions are such that p is 
near the typical value of 1.5 and is therefore too large for 
intermediate shocks to form. The magnitude of the IMF is 
increased for the second computation so that P is in the range 
for which intermediate shocks may form. All physical values 
except the IMF magnitude are the same for the two simuia- 
:ions. The solar wind speed is large enough that M K is in the 
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appropriate range tor intermediate shock formation at the 
lower value of p. The values that are the same for the two 
studies are V sw =550 km/sec. n, w =5 cm \ T sw =10 5 °K. and 
M,=10.6. The values for the high p study are B 0 =4.8 nT, 
M b =11.7, (3=1.5, and those for the low (3 study are B o =50 nT, 
M b =l.l, (3=0.014. 

The density contours for the low-beta case when an 
intermediate shock is not formed are given in Figure 3. The 
solution has been run long enough that the shock, which is 
now a fast shock, has come to an equilibrium position ahead 
of the planet At large distances from the planet the shock 
becomes very weak. 

The solution in Figure 3 is in sharp contrast to that for the 
lower value of beta shown in Figure 4. This solution is not in 
dynamic equilibrium (see below) and has evolved for a total 
of 122 Alfvdn times since the magnetic field was introduced 
into the solution. The Alfv£n time is determined from the 
Alfv6n speed computed using the above parametric values 
and the Venus radius. The bow shock has now taken on an 
almost planar shape with a depression near the Sun-Venus 
line that is concave away from the planet. This shape for the 
shock front is characteristic of one containing an intermediate 
shock in the concave portion as shown in the schematic in 
Figure 2. 

The density rise at the leading edge of the disturbance near 
the Sun-Venus line in Figure 4 divides away from the Sun- 
Venus line into one density increase due to the bow shock 



Fig. 3. Density contours for the high beta simulation dis- 
cussed in the text. The density is normalized to the constant 
initial value, and the contours are from 0.3 to 3.3 in incre- 
ments of 0.3. 



Fig. 4. Density contours for the low beta simulation discussed 
in the text (see the caption for Figure 3 ). The density con- 
tours are from 0.5 to 3.9 in increments of 0.2. The dashed 
box on the density contours locates the boundary used for the 
plot of magnetic field lines in Figure 5 . The small Hat portion 
of the density contours near the Sun-Venus iine is indicative 
of the finite grid spacing. 


and a second that extends downstream. The second rise is just 
a remnant of the original increase due to the fast shock as in 
Figure 3. The density distribution behind the planet also 
changes dramatically as the solution changes from that in 
Figure 3 with the usual fast shock solution to that in Figure 
4 with an upstream propagating shock. The larger density 
downstream of the planet in Figure 4 is due to material that 
remains there as the downstream velocity decreases when the 
shock begins propagating upstream. 

The magnetic field lines within the dashed box marked on 
Figure 4 are shown in Figure 5. The plot routine used to 
make this plot simply fits curves through local values of the 
magnetic field and arbitrarily begins and ends the curves so 
these curves should not be interpreted as being actual 
continuous field lines. The approximate location of the shock 
front, as determined by the density rise in Figure 4. is 
sketched on the field lines. Along the dashed portion of the 
shock front the shock clearly has the effect on the magnetic 
field lines expected for an intermediate shock (the component 
parallel to the shock front reverses direction). The shock is a 
fast shock along the solid part of the indicated shock front, 
and a switch-on shock separates the intermediate and fast 
shocks at the approximate location of the solid dot. 

As mentioned above, the solution in Figure 4 is not in 
dynamic equilibrium. The shock front is travelling upstream 
and is slowly decelerating. For the 24 Alfvdn times before the 
time used for Figure 4, the shock was propagating at an 
average speed of 20 km/sec, and for 73 Alfv6n times prior to 
that it was propagating at an average speed of 25 km/sec. It 
is not clear why the shock travels upstream, although it is 
probably related to the fact that the intermediate shock may 
not be as capable of reducing the solar wind momentum as a 
fast shock, and consequently the stand-off distances increases. 
The numerical box has not been made large enough to allow 
the shock to reach an equilibrium position, if it indeed 
eventually will. Such issues are naturally important to fully 
understand the physics of these solutions, although their study 
would have more relevance for the case of flow over a sphere 
than that over a cylinder. 

At large distances from the Sun-Venus line in Figure 4, 
the shock becomes a fast shock and is much weaker than it 
is closer to the Sun-Venus line where it is an intermediate 
shock. Such a weak shock at large distances is in general 
agreement with the weak shock observations reported in 
Russell and Zhang [1992]. The fast shock primarily deflects 
the magnetic field and velocity, while the intermediate shock 
compresses and slows the solar wind. 



Fig. 5. Representative magnetic field lines within the dashed 
box on Figure 4. The curve indicates the approximate location 
of the shock as determined by the leading edge of the density 
increase. 
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Discussion 

Studies ot the formation of MHD shocks near the leading 
edge ot CMEs have indicated that the type of shock formed 
depends on the plasma beta upstream of the shock and the 
magnetosonic Mach number ot the shock relative to the 
upstream flow speed [e.g„ Steinolfson and Hundhausen. 1990: 
Steinoitson. 1992]. By analogy, any dependencies on the 
parametric values found for shocks in CMEs should also exist 
in the case ot bow shocks. 

The CME studies suggest that the parametric values (beta 
and magnetosonic Mach number) during the distant Venus 
bow shock encounters reported by Russell and Zhang [1992] 
may be such that the more usual fast MHD bow shock near 
the planet is not a possible solution. Instead, we hypothesize 
(subject to two caveats discussed below) that the fast shock 
should be replaced by a shock tront containing an intermedi- 
ate shock near the Sun- Venus line and a fast shock at large 
distances from the Sun- Venus line with a switch-on shock 
separating the intermediate and fast shocks. Furthermore, the 
shock tront should be concave away from the planet over the 
portion of the front containing the intermediate shock. 

Having hypothesized that such a change in the shock types 
and configuration may occur when the solar wind conditions 
are similar to those during the observed distant bow shocks, 
we performed global MHD simulations to investigate this 
conjecture. As expected, the shock front consisted of the 
multiple shock types (switch-on, intermediate and fast) with 
the anticipated concave shape. The composite shock front 
propagated slowly upstream, and it is because of this move- 
ment of the shock front that the distant bow shock observa- 
tions can be explained. From the simulations performed to 
date, it is not clear if the shock is simply moving to a new 
( arge) stand-off distance. The shock appears to be decelerat- 
ing as it travels upstream, and the numerical box may not 
have been large enough to allow the shock to reach a 
dynamic equilibrium position. In this parameter regime it is 
evident that the bow shock location is extremely sensitive to 
relatively small changes in the solar wind conditions. Addi- 
tional simulations are needed to examine this sensitivity as 
well as to better understand the physics governing the 
upstream propagation and whether an equilibrium stand-off 
distance is always obtained. 

Although the distant shock behavior reported by Russell 
and Zhang [1992] and that in the present simulations are 
qualitatively similar, there are at least two distinctions that 
deserve comment- First, our simulations are two-dimensional 
while the solar wind interaction with a planet is clearly, in the 
global sense, a three-dimensional phenomenon. However, the 
essential physics responsible for the transition from the usual 
last shock solution to one with intermediate shocks near the 
point ot parallel propagation is not three-dimensional. It is 
certainly true that the global ramifications of this shock 
transition would differ in a three-dimensional simulation from 
that computed in a two-dimensional such as that used here. 


The important point, though, is that the necessary physics is 
included in our simulation. 

Second, at least two of the shocks reported by Russell and 
Zhang were quasi-perpendicular with angles between the 
shock normal and upstream magnetic Field of 61° and 79° 
(based on the copianarity assumption), while the upstream 
propagating bow shock in the simulations ts quasi-parallel. As 
described in the paper, local quasi-parallel propagation is 
required to make the transition to the traveling shock front 
containing an intermediate shock. Even in the more realistic 
case with an azimuthal component of the magnetic field, there 
would still be some location on the shock envelope where the 
bow shock is quasi-parallel and all of the logic and results of 
the present study would apply. The angle between the shock 
normal and the upstream magnetic field at any location 
naturally depends on local conditions. It is certainly conceiv- 
able that the necessary quasi-parallel propagation exists at one 
location and that the shocks could still be quasi-perpendicular 
at the locations reported on by Russell and Zhang. The 
present work must be viewed as basically a feasibility study, 
and more extensive simulations including the third dimension 
are needed to examine the generality and applicability of our 
results in more realistic global environments. 
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Abstract. We have developed a three-dimensional, single-fluid, MHD code which has 
successfully simulated the major qualitative features of the Venus - solar wind interac- 
tion. A bow shock forms at a height above the ionosphere close to that deduced from 
Pioneer Venus Orbiter (PVO) observations. The bow shock shows an asymmetry in 
the terminator plane qualitatively similar to results obtained from PVO data analyses. 
The magnetic field drapes over the planet and forms a two-lobed magnetotail. The 
magnetic field configuration near the planet is very similar to that recorded by PVO as 
well. We model the Venus obstacle as a hard, conducting sphere. We simulate the case 
in which the IMF is tilted with respect to the incoming solar wind velocity at approxi- 
mately the Parker spiral angle; this is the first simulation of this particular case which 
includes all spatial regions of the interaction. 


1. Introduction 

The observed features or Venus's interaction with the 
solar wind have noc yet been treated together in a single 
comprehensive model. Most models of the interaction 
have treated it as a gasdynamics process and have sim- 
plified the structure of the plasma flow in the tail re- 
gion by attaching an infinite -extension*' to the Venus 
obstacle [Spreiter and Stahara . 1980: Moore et al. . 
1991] . These approaches are a necessary first step in 
understanding the interaction. However, gasdynamics 
models have severe limitations. Thev wiil never repro- 
duce onserved phenomena such as the magnetic barrier 
behind the bow snock [Zhang ct a/., 1991] or magnet- 
ically produced asymmetries in the terminator plane 
shock configuration [e.g. Russell et al.. 1988]. Also, 
effectively extending the nightside of the Venus obsta- 
cle makes impossible the study of the near-tail region, 
where we can expect the solar wind flow to be most com- 
plex and where magnetic reconnection is most likely to 
take place. 

A model which offers improved understanding of the 
tail region of Venus i as weil as other nonmagnetized 
bodies impacted by rhe solar wind) has been devel- 
oped by Luhmaim and Schwmgenschuh [Luhmann and 
S chwi ngenschuh. 1990: Luhmann. 1990: Luhmann >t 
al.. 1991]. Here the flow outside the tail region is cal- 
culated from a gasdvnamics model like those referred 

< opyncht 1!M5 by the American Geophysical Union. 

Paoer number P5J AO'JI 7 1. 


to above, while the flow inside the tail is calculated 
by fitting a model of a comet tail inside the obstacle. 
This model includes effects of mass loading from the 
Venus ionosphere. It gives good agreement with P ho- 
bos measurements of the magnetic field around Mars 
[Luhmann et al ., 1991] at distances outside the solar 
wind bow shock. However, the above mentioned effects 
of the magnetic field on the interaction still cannot be 
studied. Also, a single comprehensive model for the 
entire interaction is still desirable. 

The first three-dimensional MHD simulation of rhe 
interaction is that of Wu [1992], YVu simulated the bow 
'hock of a hard, conducting sphere in an MHD fluid. A 
magnetic barrier was seen to form, and slight asvmme- 
rries in the bow shock caused by the interplanetary mag- 
netic field (IMF) were detected. However, the plasma J 
: i.e.. ratio of the thermal pressure to magnetic pressure) 
in this simulation was unrealistically high ( .1 zs 8). the 
rail region was not simulated, and no account was taken 
of mass loading. 

The first three-dimensionai MHD simulation of the 
V Vnus-solar wind interaction which treated all spatial 
regions of the interaction in a single model was Tanaka s 
[1993] recent simulation. Like Wu (19921. Tanaka treats 
a hard, conducting sphere in an MHD fluid. A magnetic 
barrier is formed, the magnetic field clearly drapes over 
the planet, and complex flow in the near-rail region is 
observed. 

In this paper we report the results of simulations we 
iiave made with a rhree-dimensional MHD code. Like 
Wu [1992] and Tanaka [1993]. we model the Wnus iono- 
sphere as a hard, conducting sphere and the soiar wind 
as a simrie- fluid MHD plasma. Some important physics 
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of the problem is neglected in this approximation. We 
■vili not see effects from mass loading in the rail \ Brace 
et al.. 1987], and the configurations of the shock and 
magnetic barrier can be reproduced with onlv limited 
accuracy without the ionosphere. However, this hard 
sphere model is an improvement over extending the tail 
of the obstacle, since it allows study of the near-tail re- 
gion. In the future this approach can be supplemented 
with mass loading to give a better approximation of the 
ionosphere. 

Our computational method differs greatlv from 
Tanakas [1993). He evolved the MHD equations on 
an unstructured mesh, while we use irregular spherical 
coordinates. Our time-stepping schemes are quite dif- 
ferent as well. This, in and of itself, gives our code and 
simulations value: it is desirable to have more than one 
numerical realization of the same physical phenomenon, 
since numerical studies are often fraught with hazard 
and doubt. Further, we have improved the resolution 
near the surface of the planet at the price of restricting 
the size of our computational volume. This restriction 
in volume does not significant Iv affect the results we ob- 
tain. as will be discussed below: the oniv real sacrifice 
' P ^ ave made here is that we are unable to study the far 
downstream region. We also simulate the Venus-solar 
wind interaction in the case where the IMF is tilted at 
an angle of 45° with respect to the incoming solar wind 
velocity. This case was studied by Wu (1992). but. as 
stated above, the value of J was unrealistically high and 
no attempt was made to include the tail region. This is 
the first time that all spatial regions of this asymmetric 
case have been simulated and that a realistic value of 3 
has been used. 


2. Numerical Method 
2.1. Algorithm 


_We use me modified Lax-Wendroff scheme developed 
by Rubin and Burstem (19671 to time step these equa- 
tions. A fourth-order smoother and a smoothing term 
suggested bv Laptdus (1967] are added to inhibit numer- 
ical instabilities. 

We evolve the equations in spherical coordinates with 
the # = 0 axis directed toward the Sun. Spherical coor- 
dinates allow us to pack grid points more densely nea; 
the planet s surface, where the most dramatic chang,- 
m the fluid flow occur. Near the planet's surface rh 
gnd spacing is about 0.017 R v . The grid spacing in 
creases smoothly so that about 30 grid points lie be- 
tween the planet's surface at 1 R v and a radius of 1.5 
Rv. The <y=0 axis (the Cartesian r axis) points from 
the center of Venus toward the Sun. The y axis (i.e.. 
0 = t/2. o = — /2) lies in the direction of v x xB w , vj 
and B^o being the incoming, unperturbed solar wind 
velocity and magnetic field, respectively. The x axis 
' ^ tt/ w . o 0) completes the right-handed coordinate 
system. (An alternative way of stating this is that the x 
axis is aligned so that and B-* are both contained 
in the x-z plane, and the y axis completes the right- 
handed coordinate system, i Tliis coordinate syrem wn 
by far the mosr straightforward to code and is certain! v 
the system of choice for giving a clear explanation of our 
algorithms. Unfortunately, it is not the coordinate sys- 
tem in which most observations of the Venus-solar wind 
interaction have been expressed in the literature. The 
results of our computations will be presented therefore 
in different coordinates, as will be discussed below. 

The momentum conservation (2) requires some extra 
attention. The term 

T^-(BB) 

An 

can be written 


Our simulation rode soives the three-dimensional 
nonlinear, ideal MHD commons: 

Op 

Tt + v • (pv) = .) (1) 

d(pw) _ 

~Q t 3- V - (pw) + Vp 

- — [V " BE) _ Ty-B^] =0 , 2) 

<9B 

— — V x ( v x B ) = ( t , 3 j 


-T'B TB+BV B). 

B^~ ' B is identically zero in nature. In numerical sim- 
ulations such as ours, however. V B can become large 
enough to produce sizable, nonphysical forces directed 
along the magnetic field lines. 

We have handled this difficulty by periodically sub- 
tracting the part of the magnetic field which contnbut.es 
to V - B. Tliis is accomplished as follows: the magnetic 
field can lie thought of as composed of two pieces. 

B = B f , + B-y 


0 pv 1 


p 

■ - 1 
-r- V ■ 




T'-’ )V 




v x B)j \ = 


4) 


where p is the mass density. i: , 
v is rhe thud velocity, and B is ? 


rile Thermal pressure, 
e* maimer ir held. 


where V x B = V x B„ and V B = V B,,. The held 
Brf is therefore rhe gradient of a potential <5. 

V‘$ = V B.y = V • B. 

Ue solve this equation for >f>. then find B., bv subtract- 
ing V<I> from B: 

B -+ B„ = B - VT. 

The potential T is obtained by a successive ,\or- 
. elaxation ' 1 ) R i method that we deveioned bv a< lant- 
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ing the SOR metnod of Press et ai. f 10921 to spherical 
coordinates with irregular grid spacing. 

The values or all quantities at the inflow bounaarv 
(0 < 90°) are set by input and held fixed throughout 
the time of the simulation. Quantities on the outflow 
boundary (0 > 90°) are determined by first-order inter- 
polation from the interior in such a way that djdz = i) 
for ail quantities t c representing the Cartesian axis lying 
along the 0 = 0 line ). The outer radius of the simulation 
box is 3 R\ . 

The plane containing the incoming velocity and mag- 
netic field vectors is a plane of reflective symmetry, so 
we need simulate no more than one hemisphere of the 
entire spherical volume of the system. Reflective bound- 
ary conditions hold therefore at o = 0 and o = t. 
When we are simulating flow exactly perpendicular to 
the IMF . we need simulate only a quarter sphere. In 
this case, o = 0 is still a reflective plane. At Q = t/ 2. 
symmetry boundary conditions hold for the scalar quan- 
tities. the velocities, and B ^>. while antisymmetrv con- 
ditions hold for Bo and B r . 

The Venus obstacle is modeled as a hard, conducting 
-phere sitting at the radial center of the simulation box. 
At the Venus suriace. r = 1 R\ and radiai velocities ana 
magnetic fields are set equal to zero, while all other 
quantities are determined by zero-order extrapolation 
in the radial direction. 

Our numerical method, along with most numerical 
methods, will not work unaided at the coordinates 0 = 0 
or 0 = 7T. since these coordinates are singular. We have 
found a satisfactory means of dealing with this difficulty. 
Even though this problem will present itself in almost 
any effort to model fluid or MHD physics in spherical 
coordinates, it is not extensively discussed in the lit- 
erature. We will therefore discuss in some detail our 
handling of the singular coordinates. We will discuss 
rhe 0 = 0 coordinate here: the analogy with 0 = - ran 
be drawn readily. 

Let /. j. and k label nodes in the r. 0. and o directions, 
rospecriveiv. The 7 = 1 coordinate line corresponds to 
0 = 0. For a given node on the j = 1 axis. (A l.k). the 
values ot p.p. v r . and B r are set equal to the averages 
of these quantities across the surrouding nodes (A 2. Ad. 
That is. 


; )n the axis the magnetic field and velocity lie en- 
tirely in the x-z plane which, as discussed above, is a 
symmetry plane. Therefore, on the 0 = 0 axis, there 
are no components B v or v v : B = B T x + B z z and 
v = r x x+ c-z. We have already indicated how u, and 
B z are calculated, since z = r on the 0 = 0 axis. The 
x components of these vectors A x are calculated from 
rhe corresponding vectors of the j = 2 grid points: 

T 1 V^ 1 

Ar i. i,k — ».2./02'X + A 0 t.2.iO/*X 

1=2 

1 ^ 

— ^ _ - £ Ae ,-. 2 jcos0o coso/ — A 0 i 2 .tAnot. 

1 = 2 

Now .4^ and A 0 are quite naturally defined on the 
axis as 


A» t,i ,k = A r cos < 2 * 


( 6 ) 


— — A t sin Ob. (7) 

A factor of cositf = 0) is implicit in the above formu- 
las. So. on the 0 = rr axis. 

i.i.k = -Ax cosot (8) 


•hi.U' = A t sin o^. (9) 

since cos|7r) = -1. With the components of B and v 
defined in this manner. B and v at any location along 
rhe axis are unique vectors whose components take on 
different values depending on the direction from which 
rhe axis is approached. 


2.2. Initial Conditions 


The density and pressure are initially constant 
Throughout the simulation volume. The initial velocity 
field is the flow ot an incompressible fluid over a sphere, 
directed along the 0 = 0 axis: 


<V 


— r 0 




cos 0 
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lx- 1 

r d, 2.1 1 

ik : < 

!>■ i.k 

I 'V.i.v 1 

= T 
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The k = 1 and k = A grid points are left out of the 
average, since they are boundary points and lie out- 
side of rhe proper region of computation. That is to 
sav. we wish to average over one hemisphere (or quar- 
ter sphere, in rhe case of perpendicular flow;. The co- 
ordinates f>> and O/v-i lie just inside the hemisphere 
:or quarter sphere i. while the coordinates o t and o A 
lie just outside. To include these outside points in the 
average wonid give added and uneven weighting to rhe 
aiues ar rhe o boundaries. 


= >'o 




sin 0 


l'o — 0 - 

The magnetic field is initialized in the same way. except 
that the field configuration is rotated about the y axis 
bis defined above) so that the incoming IMF is skewed 
with respect to the solar wind flow by the desired an- 
gle. With the magnetic field initialized in this manner. 
V • B is zero everywhere, discounting the effects of the 
discrete grid. The parameters we have chosen for the 
"oiar wind are summarized in Table 1. These parame- 
ters are all close to the typical solar wind parameters 
reported by Phillips rt al. [1986]. The average solar 
wind cone angle \ i.e.. t lie acute angle between rhe IMF 
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xduit; x. simulation .'solar 
Parameter 

\\ md rarametei 
Value 

-Number density 

20/cm " 

Thermal pressure 

<3.06 x 10 _i nPa 

Sound speed 

oo km/ s 

Velocity 

370 km/ . 

Sonic mach number 

6.73 

Dynamic pressure 

4.10 nPa 

Magnetic field magnitude 

10 nT 

Plasma beta ( J ) 

1.0 


md the solar wind flow direction i is given as 42.2° by 
Phillips er al.We pertorm simulations with cone angles 
of 0°. 45°. and 90°. 

3. Results 


In this section we report some of the results from our 
simulations and compare them with published analyses 
of Pioneer \ onus Orbiter (PVO) data. We give partic- 
ular attention to the configuration of the overall shock 
'tructure and to r he configuration of the magnetic field. 
- 11 though, a lew words need to be said regarding 
nomenclature. We have simulated bow shocks arising 
in IMF cone angles ot 0°. 45°. and 90°. These simula- 
tions will be referred to as the parallel. 45°. and perpen- 
dicular simulations, respectively. Also, in the following 
sections, ‘altitude refers exclusively to distance mea- 
sured from the planet's surface, while "radius'* denotes 
distance measured from the center of the planet. Fi- 
nally. up to this point we have used standard Carte- 
sian and polar coordinates ro denote spatial directions. 
From hereon, however, in order ro facilitate compar- 
isons with existing data and data analyses, we will use 
rhe B-v coordinate system of McComas ct a i. [1986]: 
rhe unperturbed solar wind Hows in the -x direction 
i corresponding to the -z direction in our code), while 
r he ct oss- Mow component ot the unperturbed magnetic 
::eld lies in rhe y direction i our code's -2 direction i. 

iie z direction lour code s — y direction i completes the 
light-handed coordinate system. 

3.1. Diagnostics 

We simulated parallel fast shocks with the three- 
dimensional (3-D) code, then compared the results with 
Emulations horn a well-tested, otten used two- 
dimensional (2-D) spherical geometry code (that is. a 
■ode written for spherical geometry and cylindrical sym- 
metry i with identical intiai conditions and spatial grid- 
ding. The results were virtually identical, as rhev must 
>e if the 3-D code is working as well as riie 2-D code. 

i. lie siiock jump conditions along the '■u agnation 
'ireamiine can be calculated analytically for rhe parai- 
cd ana perpendicular shocks. We find rfiat almost ail 
• imp conditions are satisfied in both simulations with a 
icJative error ot less than 10 1 /. The uncertainty in rhe 
ump < onditions arises trom uncertainty m determining 
fie exact location ot rhe shock and from uncertainty 
n the applicability of rhe standard shock iumo conui- 
ons. Our simulated shocks have unite wwiriis or about 


).05 /? \ ana smooth, albeit rather sudden, transitions 
between rue upstream and downstream sides. As a re- 
sult. the values of the physical quantities immediately 
downstream of the shock are often difficult to deter- 
mine exactly. The standard shock jump conditions, on 
the other hand, are derived for shock fronts separating 
two fluid regions which are in equilibrium, i.e.. two re- 
gions in which d/dt is zero for all quantities. In tlv 
system we are simulating, d/dt is zero only in the in 
flow region. Significant pressure and magnetic force* 
exist behind the shock. So even though we have rempo- 
i ai equilibrium (that is. D/Ot is zero for all quantities i, 
l/dt is certainly not zero in between the shock front 
and the planet’s surface. 

As mentioned above, a finite divergence in the mag- 
netic field gives rise to unphvsical forces directed along 
the magnetic field lines. Only if the magnitude of the 
physical force. | - Vp + B ■ VB/4 tt - VB 2 /Sn\. is much 
larger than the magnitude of the purely numerical force. 
!BV • B/4tt|. can we believe the results of our simula- 
tions. We have inspected the ratio of these two forces 
-n both the perpendicular and 45° simulations. We find 
'hat almost everywhere in the simulation volume rii 
numerical force is less than V/ of the physical force. 
The exceptions to this rule occur in isolated volumes 
within the numerical shock front and result from the 
wav our numerical procedure handles discontinuities. 

We have examined how the location of the computa- 
nonal outer boundary affects the simulation results by 
performing simulations of the parallel shock with the 
outer boundary set at 3 and at 10 R\ . .V comparison 
of pressure contours in the region extending from the 
surface out to 3 R\ is shown in Figure i. Figure 1. top. 
ffiows the pressure distribution with the outer computa- 
tional radius set at 3 Rv * and Figure i. bottom, shows 
rhe distribution when the outer computational radius is 
at 10 R\ . Differences are apparent but small. Likewise, 
dl other quant ires show only small differences which 
\re in no wav qualitative. On rhe basis of rhis study. 
vp have < uosen to place rhe outer radial boundary tor 
ail the present simulations at 3 R\ in order ro obtain 
better numerical resolution near the planer’s surface. 


3.2. Global Behavior 

Figure 2 shows field lines from the final configuration 
ot the 45° simulation. The line traces were begun at 
points on the outer computational boundary. 0.1 R\ 
above tiie ecliptic plane. The held lines piie up at rhe 
mock front, then slip over the planet. W irhin the siiock 
rhev are advecred downstream less quickly than outside 
rue shock. Bv rhe time they have escaped the vicintv or 
rhe planer s mriace on rhe downstream side, rhev are 
a retched ana hem so that they form a rwo-lobed mau- 
uetotaii behind the pianet. The magnetic field points 
largely sunward in one lobe (the top iobe of Figure 1. in 
rhis case) and largely anrisunwani in the other. The two 
lobes are separated bv a distinct current sheet [McCo- 

;nas ft ai. 19361. This behavior qualitatively matches 
'he behavior predicted bv A If von (19571 and observed bv 
• 'Corn as - r ; 19861. ?Wre rtiat iu me bow siiock tin* 
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contour to oe gee :cntour interval of i iin pto.si- 2i.iaB 

Figure 1 . Pressure contours of the solar wind impacting Venus when incoming flow and inter- 
planetary magnetic field (IMF) are aligned along the sunward axis. The top half sho\vs the results 
of a simulation performed with the outer computational boundary set to 3 Ry: the bottom half 
show's the results of a simulation with the outer computational boundary at 10 Ry. 



Figure 2. Magnetic field lines of the solar wind plasma 
impacting \ onus. The solar wind flow's in from the left 
and out the right. The incoming magnetic field is rilted 
at an angle of 4 o° with respect to the incoming velocity. 
Note that the srrerched magnetic field lines form a rwo- 
iobed magnetoraii behind the planer. 


magnetic field line piieup is greatest where the field is 
more closely perpendicular to the shock normal. This 
is expected, since, all other things being equal, as a 
shock becomes more perpendicular, the increase in the 
magnitude of the magnetic field becomes larger. The 
perpendicular simulation shows all of the above global 
characteristics, but, of course, it is symmetric with re- 
spect to the x-r plane ti.e.. the .r-r plane of the B- 
r coordinate system). (A note of caution: magnetic 
field magnitudes cannot be reliably inferred from the 
field line densities in Figure 2. First, perspective effects 
make some lines iook closer than thev actually are. Sec- 
ond. the starting points of the lines, particularly in the 
rail region, were not chosen to systematically indicate 
magnetic field magnitude but to clearly show the global 
behavior of the field lines, i 

It can be seen in Figure 2 that the field lines of the 
unperturbed magnetic field are not perfectly straight. 
This is because of the intiai conditions discussed in the 
-ection 3.1. The initial upstream field eouid not be 
• hosen to be perfectly uniform because we needed to 
-atisfv V ■ B = ti in the presence of the conducting 
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cross the surface of the planet. In the simulation that produced these results the IMF was 
perpendicular to the incoming velocity. Each of the five elements was given an initial SZA of 
about 0 and initial polar angle positions of 0°. 18°. 36°. 54°. and 90° measured from the ecliptic 

K, CL i C / * Un *g« t forces manifest themselves, in that fluid elements moving in the noon- 
mdnight plane (o - 90 ) are accelerated more quickly than fluid elements moving in the ecliptic 


planet s surface. Therefore the initial magnetic fields 
are not exactly parallel to each other. Since the inflow 
boundary conditions are held fixed throughout the run 
of the simulation, the incoming upstream magnetic field 
are slightly curved. Setting the outer computational 
ooundarv farther out wouid ameliorate this effect. The 
slight departure from uniformity in the upstream mag- 
netic field seems to have no significant effect on the 
plasma flow. 

Magnetic field slingshot forces manifest themselves in 
the velocity distribution of the plasma near the planet's 
surface. Figure 3 shows, as a function of time, the solar 
zenith angle of fluid elements in the perpendicular simu- 
lation moving across the surface of the planet. Here the 
initial positions of the fluid elements were chosen just 
above the planet's surface, near the subsolar point. The 
plasma moving across the pole (© = 90°) is accelerated 
more quickly and reaches the far side of the planet faster 
than the plasma moving in the ecliptic plane. This hap- 
pens because the perpendicular field lines become "hung 
up” in the bow shock, then slip quickly over the pole, 
dragging solar wind plasma with them. 

3.3. Shock Configuration 

We find that the standoff distance ar rhe subsolar 
point varies with the cone angle, but we do not see a 


monotonic trend. If we locate the shocks by the outer 
edges of each shock layer, the standoff distance of the 
parallel shock is 0.28 R Vm the standoff distance of the 
45° shock is 0.22 J?x*. and the standoff distance of the 
perpendicular shock is 0.26 Rv. In absolute terms these 
numbers translate into 1700. 1340. and 1580 km. re- 
spectively. The 45° simulation is the one most repre- 
sentative of the typical solar wind conditions at Venus 
[Phillips et al . . 1986]. Zhang et al. [1990] find from 
PVO data that the subsolar standoff distance varies 
from 1600 km at solar minimum to 2200 km at solar 
maximum. The difference between our computed re- 
sults and the observations can be explained by the lack 
of a model of V enus s ionosphere in our simulations. If 
we try to account for the ionosphere by simply adding 
300 km to the effective radius of Venus (approximately 
the altitude found by Zhang et al. [1990]). we can ex- 
pect the computational standoff distances to scale pro- 
portionately. since the Venus radius is the only inherent 
^caie length in the system. We then find standoff dis- 
tances above the ionosphere of 1780. 1400, and 1650 km. 
To find the altitudes of the shocks above the planer it- 
self. we add in the 300 km thickness of the ionosphere. 
This gives altitudes of 2080. 1700. and 1950 km. These 
numbers are now well within the range found by Zhang 
et al. Including the ionosphere, then, is essential for de- 
termining the bow shock location. 
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Figures 4-8 show distributions of pressure and mag- 
netic field magnitude from the parallel. 45° . and per- 
pendicular simulations. The outer computational bound- 
ary was 3 R v% but here we show only the region con- 
tained within 1.5 Rv in order to bring out more details. 

The pressure and magnetic field magnitude distribu- 
tions of the parallel case are shown in Figure 4. Since 
the parallel configuration is cylindricaily symmetric, the 
distributions shown here represent the distributions in 
ail polar angle planes. The decrease in the magnetic 
field strength downstream of the shock, around the 
subsolar line, is the result of MHD physics together 
with the boundary conditions imposed on the mag- 
netic field. That is. a parallel shock has no increase 
in magnetic field, so the field strength remains fairly 
constant through the shock. Also, the field decreases 
near the pianet ? s surface because the planet is conduct- 
ing: there can be no radial magnetic field component 
at the planet s surface, but near the subsolar line, the 
magnetic field is almost purely radial. Therefore the 
magnetic field strength drops to zero in this region. 

The distributions shown for the 45° (Figures 5 and 
>i and perpendicular (Figures 7 and 8) cases represent 
only the ecliptic and noon-midnight pianes. The mag- 
netic field strength rises monotonically to a peak at the 
planet s surface, while the pressure reaches a maximum 
just downstream of the shock, then drops to a lower 
value at the planet's surface. This behavior indicates 
the formation of a magnetic barrier, which we discuss 
in more detail below. 

The solution for the 45° simulation does not pos- 
sess reflective symmetry about the x axis l i.e., the axis 



(a) 



Figure 4. iai Pressure and (b) magnetic held magni- 
tude distributions of the solar wind. The maenetic field 
and the incoming vriocitv are aliened alone the sunward 
axis. 



(a) 


(b) 

Figure 5. Pressure distributions of the solar wind with 
(a) the ecliptic plane and (b) the noon-midnight plane 
shown. The IMF is tilted by 45° with respect to the 
incoming velocity. The IMF orientation is such that the 
shock is quasi-perpendicular in the top half and quasi- 
parallel in the bottom half of Figure 5a. 


pointed toward the Sun or the 0 = 0 axis in spherical co- 
ordinates), a symmetry which the parallel and perpen- 
dicular solutions possess by definition. For instance, the 
part of the shock front where the shock normal is quasi- 
perpendicuiar to the magnetic field stands about 6% 
farther out from the planet's surface than does the part 
of the front where the shock normal is quasi- parallel. 
This is to be expected: the magnetosonic speed along 
the shock normal in the quasi- perpendicular region is 
higher than the speed in the quasi- parallel region. Con- 
sequently. information in the quasi- perpendicular region 
can propagate farther away from the planet's surface 
before its progress is impeded by the onrushing solar 
wind. This is aiso consistent with ITVs [1992] MHD 
simulations of a plasma with J=8. 

The distributions of the plasma quantities within the 
bow shock are also asymmetric for the 45° case. The 
magnetic pressure has a slightly higher buildup in the 
(inasi-perpendicuiar region than in the quasi- parallel 
region. This is consistent with the general expecta- 
tion that perpendicular shocks have a larger buildup 
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(a) 



(b) 

Figure 6. Magnetic field magnitude distributions of 
the solar wind with (a) the ecliptic plane and (b) the 
noon-midnight plane shown. The IMF is tilted bv 45° 
with respect to the incoming velocity. The IMF orienta- 
tion is such that the shock is quasi-perpendicular in the 

top half and quasi-parallel in the bottom half of Figure 
6a. ° 


of magnetic field than parallel shocks. The thermal 
pressure and density maxima shift away from the sub- 
solar point toward the quasi-parallel side of the shock. 
The slightly higher buildup of magnetic pressure in the 
quasi-perpendicular region might explain this, since it 
would be expected that this pressure would divert some 
of the plasma flow into the quasi- parallel region. 

In the terminator plane we find average distances 
from the planet center to the bow shock of 1.76. 1.92. 
and 1.95 Ry. for the parallel. 45°. and perpendicular 
shocks, respectively. The terminator locations of the 
45° and perpendicular shocks are shown in Figure 9. 
The locations have been determined by the rise in en- 
tropy of the shocked solar wind. Regarding the latter 
two shocks, there is an uncertainty in their positions of 
the order of 0.1 Ry because, as mentioned above, the 
simulated shocks are not true discontinuities. The par- 
allel shock, on the other hand, is almost as well defined 
at the terminator as it is at the nose. Russell et al. 
[1988] find that the average terminator distance varies 
between 2.14 Ry at solar minimum and 2. 10 R v at solar 


maximum. The discrepancy here is. again, attributable 
to the lack of the Venus ionosphere in our simulations. 

We find that the average shock altitude at the termi- 
nator increases as the cone angie increases. This agrees 
with Zhang et al. [1990], who find that the altitude 
is typically 1630 km for cone angles less than 45° and 
1840 km for cone angles greater than 45°. The relative 
difference is about 12%. The relative difference betwee • 
our perpendicular and parallel altitudes is around 25'. 
This is not necessarily inconsistent with Zhang et a, 
As they state, if they could have had more data at low 
cone angles, their relative difference would have been 
much larger. On the other hand, we have compared the 
two most extreme cases possible. Comparisons of shock 
altitudes corresponding to cone angles of 30° and 60° 
would be more appropriate. 

We find that the shock boundary in the terminator 
plane is asymmetric in a manner qualitatively consistent 
with the findings of Russell et al. [1988], They find that 
when the IMF is nearly perpendicular to the solar wind 
flow, the distance from the planer center to the shock in 
the noon-midnight plane is about 10'/ greater than the 
distance along the equatorial plane. Our perpendicu- 



(a) 



Figure 7. Pressure distributions of the soiar wind with 
a) the ecliptic piane and (b) the noon-midnight plane 
-mown. The IMF is perpendicular to the incoming solar 
ind. 
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(a) 



(b) 

Figure 8. Magnetic field magnitude distributions of 
the solar wind with (a) the ecliptic plane and (b) the 
noon-midnight plane shown. The IMF is perpendicular 
to the incoming solar wind. 


iar simulation gives a smaller asvmmetrv of about 5%.. 
which is consistent with Tanaka's [19931 recent MHD 
simulation result. Tliis asymmetrv can be seen, with 
some effort, in Figure 9. It has been suggested bv Rus- 
sell et ah. among others, that the asymmetry is caused 
by the difference in the polar and equatorial magne- 
tosonic speeds. Given our code's radial resolution at the 
terminator shock location (A r « 0.03). our code should 
be able to differentiate between a 5% and 10% asymme- 
try. In light of our agreement with Tanaka and in view 
of our radial resolution, it seems that the magnetosonic 
speed accounts for some but not ail of the asymmetrv. 
The rest of the asymmetry arises from non-MHD pro- 
cesses. perhaps pole-equator asymmetries in ion pickup, 
as mentioned bv Russeil et al. 

The ratios of the terminator shock radii to the sub- 
solar or nose radii are given by 1.44. 1.57. and 1.35 for 
our three simulations. Zhang 1 1 ai (1990) find a typical 
value of 1.66. 

3.4. Magnetic Fields in the Shock 

The most striking feature of the behavior of the mag- 
netic field is the dear draping of the Held as it is ad- 


vected around the planet las shown m Figure 2 and 
discussed above). The ecliptic plane magnetic field con- 
figuration of the 45° simulation is shown in Figure 10. 
The solar wind is flowing in from the top. Regarding 
the magnetotail, the magnetic field of one lobe is fairly 
uniform in direction and has a magnitude of the or- 
der of 10 nT. In the other lobe the magnetic field lines 
are strongly curved and the magnetic field strength is 
0.1- 1.0 nT. The results shown in Figure 10 are very 
similar to the average magnetic field distribution found 
by Phillips ft. ai [1986). (See. particularly. Phillips 
et al. [1986, Figure 5d]. Keep in mind that because 
of the orientation of our simulation magnetic field, our 
results must be reflected across the x~z plane before a 
direct comparison can be made with PYO data. This, 
of course, entails no change in the physics involved 
here.) Two points deserve particular mention in this 
regard. First, at low altitudes the field tends to encir- 
cle the planet. Second, the separating current sheet is 
not aligned with the plasma flow. Rather, it is angled 
toward the lobe of weaker magnetic field. Also regard- 
ing the current sheet, we find that it is displaced from 
the center of the tail by about 0.5 R\ . Me Comas et ai 
[1986) find a similar displacement 8-12 R\- downstream. 

The magnetic field magnitude from the 45° simula- 
tion is shown in Figure 11. The half of the ecliptic 



Perpendicular simulation 

45° S im U iation 

Figure 9. Cross sections of the shock locations in the 
terminator plane. The view is from upstream, looking 
downstream. The solid line is the cross section pro- 
duced by the simulation where the IMF was perpendic- 
ular to the solar wind velocity. Note the slight asym- 
metry; the distance to the shock is somewhat longer 
in the noon-midnight plane than in the ecliptic plane. 
The dashed line is from the simulation where the IMF 
was tilted at 45° with respect to the velocity: the left 
-iae of this cross section is where the shock was most 
perpendicular. 
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Figure 10. Magnetic field in the ecliptic plane. The IMF is tilted by 45° with respect to the 
incoming velocity. The solar wind plasma is flowing in from the top. The vectors in the top right 
comer represent magnitudes of 10 nT. 


plane in which the shock is quasi- parallel is shown in 
Figure 11a, the noon-midnight plane in Figure lib. and 
the half of the ecliptic plane where the shock is quasi- 
perpendicular in Figure 11c. At solar zenith angles sur- 
rounding the terminator we find magnetic field magni- 
tudes similar to those found bv Phillips * t ai [1986. 
Figures 3a and 4a]. Differences near the planet's surface 
may be caused by the ionopause slowing the flow of the 
plasma and magnetic field [Me Gary, 19931. 

Figure 12 shows the values of the sunward B x and 
cross-flow By components of the magnetic field of the 
perpendicular simulation "measured along a circular 
orbit of radius 1.4 Ry. The values of the sunward com- 
ponent are very similar in magnitude to those recorded 
b y PVO at similar altitudes and in similar solar wind 
conditions [Luhmann et ai . 1991: Dubinin et a/.. 1991]. 

The cross-flow component is another matter, how- 
ever. At the current sheet in the near-tail region we 
find that the cross-flow component of the magnetic field 
is quite small. This is consistent with Tanakas [1993] 
MHD simulations. However, this is very different from 
the results of the analysis of PVO data by Luhmann et 
ai [1991]. Their measurements show that the cross- 
flow component reaches a sizable maximum irs 10 nT) 
at the current sheet. Although this maximum is not 
always a clear feature in observations reported by other 
authors, the cross-flow magnetic field in this region is. 


in fact, consistently higher than our result [ Luhmann et 
ai, 1991: Brace et a/.. 1987: Dubinin et ai . 1991]. It 
would appear that processes apart from MHD phenom- 
ena are at work maintaining the cross-flow field strength 
in this region. For instance, this could be an indication 
of the importance of mass loading in rhe tail region. 

Inspection of Figure 12 will show not only that the 
cross-flow component of the magnetic field becomes 
small, but also that it actually changes sign as well. 
This sign change indicates field reversal in the cross-tail 
component. This is consistent with magnetic reconnec- 
tion. In fact, the field line traces of Figure 13 show 
that a region of reconnected magnetic field lies directly 
behind the planet. The field lines shown in Figure 13 
were traced from 20 starting points lying in the tail re- 
gion of the noon- midnight plane about 0.1 R v above the 
ecliptic plane. The reconnected fields are stretched out 
into the tail in a roughly conical volume with a height of 
about 1 Ry and a base radius of about 0.3 R v . It would 
not be possible to produce a coherent figure that gave a 
representative sample of all the reconnected field lines 
in the simulation. However, we have found that recon- 
nected field lines cover the surface of the planet, and 
rhat they cover the the nightside much more densely 
than the davside. 

The lower third of the reconnection region is within 
Tie altitude ot rhe niehtside ionosphere. So. if this re- 
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Figure 11. Contours of magnetic field magnitude 
showing (a) the half of the ecliptic plane where the 
shock is quasi- parallel, (b) the noon-midnight plane, 
and (c) the half of the ecliptic plane where the shock is 
quasi-perpendicular. Labels show magnitude of the lo- 
cal magnetic held relative to the magnitude of the IMF. 
The IMF magnitude is 10 nT. The IMF is tilted at 45° 
with respect to the incoming flow velocity. 


gion does indeed exist in the near tail of Venus, ir is cer- 
tainly located at a higher altitude than we have found 
here <uid is probablv much more complex in its con- 
figuration. Nevertheless, our finding that reconnection 
takes place over a very limited region of the magnetotaii 
is consistent with McComas t:t al.'s [1986] result that 
8 - 12 Rv downstream the cross-flow held is reversed 
oniv about 5/ of the time. The resistivity producing 



Figure 12. Sunward or x component (solid line) 
and cross-flow or y component (dashed line) of mag- 
netic field in a circular orbit in the ecliptic plane 0.4 
R v above the Venus surface. The cross-flow compo- 
nent becomes slightly negative in the tail region (orbit 
anglew 180°). indicating field reversal and magnetic re- 
connection. The IMF is perpendicular to the flow ve- 
locity. 


our reconnection originates in numerical dissipation and 
is not a part of our algorithm. We have found recon- 
nection only in the perpendicular simulation. If there is 
reconnection in the 45° simulation, it is occuring only in 



Figure 13. Reconnected field lines in the tail region. 
The IMF is perpendicular to the incoming flow velocity. 
The reconnected lines are pulled away from the planet 
into a roughly conical volume of height 1 Rv and base 
radius of 0.3 R\ . 







21.656 


.-DIMENSIONAL MHD SIMULATIONS 


CABLE AND STEINOLFSON: THREE 


/ 

/ 

/ 



Figure 14. Field lines showing spiraling behavior in 
rhe rail region. These lines pass rhrough rhe noon- 
midnight plane ar verv low altitudes above rhe planer 
on the nightside. The IMF is perpendicular ro rhe in- 
•■omuiff rlow velocity. Some reconnected field lines are 
mown as weil. 


cl re ?i°n sma d enough to have escaped our determined 
searching. 

We find complex behavior in the magnetic field lines 
behind rhe planet in both the 43° and perpendicular 
cases. The most sunward (and antisunwardi directed 
lines spiral around one another as they head out from 
the planet. Some of these field lines from rhe perpen- 
dicular simulation are shown in Figure 14 alone with 



Altitude < Rv) 


°me. reconnect en heid lines, i: -houid he noted rl,.„ 
-mce our simulation box extends ro <>mv 3 By. ;hese 
lines are not "anchored" in the solar wind plasma out- 
side the shock. Therefore extending the outer boundary 
of the simulation may alter this behavior. 

Looking ar rhe nose of the shock, we see rhe formation 
ot a definite magnetic harrier ar the subsolar poim 
both the 43“ and perpendicular simulations. Figure 
shows rhe values of the magnetic and thermal pressn: 
relative ro rhe solar wind rhermal pressure, along : 
’tibsolar hue in both simulations. Mb magnetic i,.... 
rier is formed in the parallel simuiaron. as is consistent 
with .MHD magnetic flux conservation and our bound- 
ary conditions. See the discussion of Figure 4. i The 
rhermal pressure of rhe 45° simulation peaks ar about 
0.08 R i above the planet's surface, reaching a value of 
about 4.12 nPa. The magnetic field rises monotonicallv 
to a peak value of only 2.85 nPa. In rhe perpendicular 
emulation rhe pressure peaks ar abour U.l B\ . reaching 
<i value of about 3.70 nPa. Ar the planer's surface the 
pressure drops to abour 3.12 nPa. The magnetic pres, 
-ure. on riie orher iiand. rises monoronicaiiv ro a van: 

>1 3.88 nPa ar rhe surrace. 1 he niauuenr nressiire in : : 
ppi penciicuicii case makes up a much larger percental, 
of the total pressure than in the 45° case. As indicated 
in Table 1. the solar wind dynamic pressure is 4.10 nPa. 

Measuring the thickness of rhe magnetic barrier is im- 
piecise not onlv because our shocks have finite widths, 
but also because rhe definition of rhe magnetic barrier 
is somewhat arbitrary. If we define rhe magnetic bat- 
her location by Zhang , t. al.' s [1991] criterion as rhe 
point where r lie magnetic pressure is <-(|ual ro half of 



Figure lo. r dermal and magnetic pressures plotted along the subsolar 
simulation and lb) perpendicular simulation. The shaded regions indicate 
.avers m earn cn.se. In both simulations, rne soiar wind uvnamic on-ssura 
•inoerriiroefi rr.ermai pressure, has a vmue ot 07.7. 


line for rhe sn. 45' J 
Mic shock boundary 
wiu>n sraied i v mo 
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Hie upstream uvnamic pressure, we ger thicknesses ot 
0.1)5 aim u.l i?v . These numbers rransiate into 300 and 
600 km. Zhang >.t, ai. [1991] find that the magnetic 
harrier has a median altitude of about 1000 km when 
the upstream dynamic pressure is around 4.G nPa (our 
present dynamic pressure). Recalling the value of 300 
km given bv Zhang * t al. [1990] for rhe typical iono- 
sphere thickness. this gives an observed magnetic bar- 
rier thickness of 700 km. If we try to account for the 
ionosphere here in the same way as we did when dis- 
cussing the bow shock location in section 2.2. we find 
magnetic barrier thicknesses of 330 and 630 km from the 
45° and perpendicular simulations, respectively. These 
numbers fall short of Zhang at ai' s [1991] value. How- 
ever. it must be remembered that the ionosphere is not 
a perfect sphere but is a rather oblate obstacle: its ra- 
dius or curvature at the subsolar point is significantly 
larger than 1 R\ . Could we account for the shape of rhe 
ionosphere in a simple manner, the thicknesses of our 
simulated magnetic barriers would increase even more, 
perhaps agreeing very well with Zhang ct al.‘s value. 

4. Discussion 

Our three-dimensional MHD code lias successfully 
simulated rhe observed major qualitative features of the 
interaction between Venus and the solar wind. Our sim- 
ulations show a how shock at slightly lower altitude 
than observed bv PYO. a magnetic barrier beneath the 
bow shock, asymmetry in the terminator cross section 
of the shock, and a two-iobed inagnetotail formed by 
draping of the held lines over the planet s surface. All 
quantities given by the simulations (except the near- 
tail magnetic field) are similar in magnitude to those 
determined from PYO data. 

We have also shown some limits of single-fluid MHD 
as a model of rhe Venus-solar wind interaction. Most 
obviousiv. rlie lack of an ionosphere changes the sc;ile 
>t nmnv nr rhe resuits such as rhe bow shock srandoff 
distance and rhe magnetic barrier rliickness. Including 
rhe ionosphere "by hand." he., adding 300 km to the 
effective radius of Venus, brings the standoff distances 
into the range of observations. However, the magnetic 
barrier thickness remains too small. The lOVf asymme- 
try of the shock in the terminator plane is only par- 
tially accounted for by single-fluid MHD. Some other 
process must be found to explain a missing 5/. Fi- 
nally. single-fluid MHD does not give correct results for 
the qualitative behavior of the cross-flow magnetic field 
component in die near-rail region. In this region we 
find cross- riow magnetic Helds which are much smaller 
than those typically observed. 

i hir simulations have produced several results which 
^eem to have nor been srudied in depth heretofore. In 
die as'vmmotne LYsirnulation the peak in the magnetic 
held buildup i-* muted from rhe subsoiar point toward 
the quasi-perpendicular region of the shock. In con- 
trast. rhe peaks in thermal pressure and density art? 
-bitted toward die uuasi-parailel region. This is ro be 
sneered. onasi-tcTprnuiniiar shocks have nigiier 


buildup qi - magnetic rieid. and such a buiidun can be ex- 
pected to divert piasma riow toward rhe quasi- parallel 
region. In the near-tail region our perpendicular sim- 
ulation has established the possibility of magnetic re- 
connection. although the magnetic field configuration 
in the actual near- rail region of Venus is certain to be 
more complex than what we have seen in our simula- 
tions. Lastly, in the perpendicular and 45° simulations 
we have seen a new class of field lines which exhibit 
a spiraling behavior as they are traced away from the 
planet s surface. This spiraling behavior seems to he 
intermediate in some sense to reconnection on the one 
hand and to smooth motion of the field lines along with 
the plasma flow on the other. 

In future research we will model the ionopause and 
add mass loading terms to the MHD evolution equa- 
tions. On the basis of the outcomes of our present sim- 
ulations. we have good reason to hope that future re- 
sults will show improved quantitative agreement with 
the PYO data. 
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Numerical simulations of mass loading in the solar 
wind interaction with Venus 

K. Murawski 1 and R. S. Steinolfson 
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Abstract, Numerical simulations are performed in the framework of nonlinear two 
- dimensional magnetohydrodynamics to investigate the influence of mass loading 
on the solar wind interaction with Venus. The principal physical features of the 
interaction of the solar wind with the atmosphere of Venus are presented. The 
formation of the bow shock, the magnetic barrier, and the magnetotail are some 
typical features of the interaction. The deceleration of the solar wind due to the 
mass loading near Venus is an additional feature. The effect of the mass loading is 
to push the shock farther outward from the planet. The influence of different values 
of the magnetic field strength on plasma evolution is considered. 


Introduction 

Data from the Pioneer Venus Orbiter (PVO), as well 
as data from Venera 9 and 10 and other sources, have 
contributed to a growing understanding of the inter- 
action of the solar wind with an unmagnetized planet 
such as Venus. The general picture derived from the 
available data is that, for typical solar wind conditions, 
an electrically conducting ionosphere deflects the super- 
sonic and superalfvenic solar wind around the planet. A 
bow shock forms upstream of the ionosphere and serves 
to slow and also assists in deflecting the solar wind. 
The size of the bow shock depends on the sonic and 
alfvenic Mach numbers and on the shape and the size 
of the effective obstacle, whereas asymmetries in the 
shock shape are determined by the direction of the in- 
terplanetary magnetic field (IMF) [Khurana and Kivel- 
son, 1 994 ] . Zhang et al. [1990] have proven that at 
Venus the effective obstacle undergoes a systematic size 
change of 0. 13 R^, (where R 0 is the radius of Venus) dur- 
ing a solar cycle. The boundary separating the shocked 
and magnetized solar wind plasma from the thermal 
ionosphere plasma is referred to as the ionopause, and 
the region between the bow shock and the ionosphere is 
referred to as the magnetosheath. The inner portion of 
the magnetosheath contains a region of enhanced mag- 
netic pressure referred to as the magnetic barrier [e.g., 

Luhmann, 1986]. There are physical processes involved 
in these relatively large-scale interaction regions as well 
as in the magnetic tail that the present work has appli- 
cation to. 
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Zhang et al. [1991] have used PVO data to show 
that the convected field gasdynamic model of Spreiter 
et al. [1966] predicts the correct bow shock location 
if the magnetic barrier is taken as the obstacle around 
which the plasma is deflected. They also found that the 
magnetic barrier is strongest and thinnest (about 200 
km thick) at the subsolar point and becomes weaker 
and thicker at the terminator. The magnetic barrier is 
responsible for transferring most of the solar wind dy- 
namic pressure to the ionosphere through the enhanced 
magnetic pressure. 

Observations at Venus have revealed that a magnetic 
tail forms as a consequence of the solar wind interac- 
tion with the Venusian ionosphere. The portions of IMF 
lines passing near the surface of Venus are slowed due 
to the interaction with newly ionized atmospheric neu- 
trals. The two ends of the magnetic ropes continue to be 
pulled downstream by the solar wind [e.g., Spreiter and 
Stahara f 1980]. Slavin et al. [1989] have estimated that 
the magnetotail extends to between 50 to 150 from 
the center of the planet for a typical solar wind flow 
with the Alfven speed V A = 60 km/'s and the Aifven 
Mach number = 7.2. 

A study by Nagy et at . [1981] showed that Venus 

has a dayside neutral exosphere dominated by oxygen 
at altitudes above — 400 km calculated from the plane- 
tary surface. The exospheric oxygen together with other 
particles present can be ionized by solar radiation and 
charge exchange. Newly ionized particles are eiectro- 
magnetically coupled to the shocked solar plasma of the 
magnetosheath. As a consequence, near the ionopause 
the plasma becomes mass loaded due to interactions 
with ions, particularly 0+ [e.g., Luhmann et a/., 19911. 
Model computations by Gombosi et al.. [1980! using in- 
put from PVO ionosheath and neutral atmosphere pro- 
files indicate that approximately 10 % by number of the 
solar wind protons may undergo charge exchange with 
neutrals at ionopause altitudes. As a result of the mass 
loading (ML), the plasma is slowed and compressed. 
When the plasma then flows around the planet into the 


2547 



2548 


MURAWSKI AND STEINOLFSON: SIMULATIONS OF MASS LOADING 


magnetosheath, momentum conservation requires an in- 
crease in flow speed into the magnetotail. 

It has been noted by Luhmann [1986] that the above 
general behavior tends to prevail providing the thermal 
pressure in the upper ionosphere is sufficiently larger 
than the solar wind pressure and the pressure balance 
occurs where the ionospheric plasma is coilisioniess. As 
the solar wind pressure increases, relative to that of 
the ionospheric plasma, the magnetic barrier eventu- 
ally breaks down and the height of the dayside iono- 
sphere decreases. When the ionopause is too close 
to the planet, the transterminator transport of iono- 
spheric plasma ceases, thereby leading to a nightside 
phenomenon known as the disappearing of the iono- 
sphere f Cravens et a/., 1982], 

Cloutier et al. [1987] have noted that ML processes 
at Venus should be asymmetric since the solar wind 
electric field depends on the relative directions of the 
flow velocity and the IMF. Phillips et al [1987] have 
claimed that the flow and field configuration of the mag- 
netosheath plasma, together with the large gyroradius 
of the pickup ions, cause ML to occur preferentially on 
one side of the magnetosheath. More efficient pickup 
of newly created ions should occur over the hemisphere 
that produces an electric field directed outward, away 
from the ionosphere. Luhmann et al. [1991] and Phillips 
et al [1986] have presented observations indicating an 
increase in the magnetic field intensity and current over 
the hemisphere where the electric field is directed out- 
ward. The enhanced solar UV during solar maximum 
increases the ionization rate of neutrals, which moves 
the bow shock out from the planet both in the equa- 
torial plane and at the terminator [ Zhang et al. } 1990]. 
Similarly, Alexander and Russell [1985] have shown that 
the position of the bow shock at the terminator de- 
pends on solar activity. This dependence may be indi- 
rect evidence of ML since, when the solar EUV is high, 
more mass may be added to the shocked solar wind, 
thereby forcing the bow shock to recede from the planet. 
Linker et al [1989] have revealed that ML can increase 
or decrease the plasma temperature, depending on the 
value of the sonic Mach number A/ s . When 7 = 5/3, 
Ms > v'9/5 is required for heating to occur. Numerical 
simulations performed show that the effects of ML on 
plasma temperature and velocity are most pronounced 
in & wake region. However, this theory was developed 
for the case of Io’s atmosphere, which is described by 
MHD equations with the source terms in the mass con- 
tinuity, Euler, and energy equations. In the case of Io, 
the magnetic field lines are in the north-south direction 
and are perpendicular to the horizontal flow. 

Recent computer simulations of the three-dimensional 
global interaction between the solar wind and unmag- 
netized bodies have proven to be extremely useful for 
improving our understanding of the associated large- 
scale processes [e.g., Wu, 1992; Tanaka , 1993; Cable 
and Steinolfson , 1995; Gombosx et al 1994; McGary 
and Pontius, 1994). Venus is of particular interest since 
the planet does not have a significant intrinsic mag- 
netic field, and the interaction of the solar wind with 
the Venusian ionosphere involves fundamentally differ- 


ent physical processes than occur at magnetized Earth. 
In addition, a large quantity of relevant observations arc- 
now available for Venus. When interpreted in associa- 
tion with the simulated results, this data set provides a 
test of the physical processes included in the model. 

Single-fluid MHD simulations of the global interac- 
tion of the solar wind with Venus without considera- 
tion of the ionosphere have been performed by several 
investigators. Wu [1992] limited his study to the day- 
side. Tanaka [1993] included the nightside as well but 
only considered magnetic field orientations parallel and 
perpendicular to the solar wind flow and was in a pa- 
rameter regime not directly applicable to average solar - 
wind conditions at Venus. Cable and Steinolfson f L 995 i 
carried out simulations for typical observed solar wind 
conditions at Venus and for an IMF orientation near 
the Parker spiral. McGary and Pontius [19941 studied 
the effects of ML, but they considered only the dayside. 
Moreover, they discussed the case of a cylinder with 
d/dz = 0 and B along the cylinder axis (z - axis) ana 
the perpendicular flow, namely, V J_ B. 

Our purpose is to extend the model of McGarv ana 
Pontius to study the nightside and to discuss various 
strengths of the magnetic field. We also considered a 
more realistic spherical model. 

The paper is organized as follows. The physical model 
used in the present study is described in the next sec- 
tion. Section 3 describes the numerical model which is 
used in the present studies. We present and discuss our 
detailed results in section 4. The paper concludes with 
a short summary. 


Physical Model 

We solve the following form of the compressible MHD 
equations as an initial value problem: 

dp 

+ V * (pV) = qo exp [-(/i - ho}/ Ho], i j 


3(pV) 

+ v ■ f^ v ) v l = 


-Vp+ — (V X B) X B. (2) 

4 7T 


0B 

dt 


= V x (V x B), 


(3) 


d.pV* P B\ r P V 2 y 

»<— + ^T + i7) + v «V-TTT« v 


-Til* * (V x B)l) = ,). 


Here p is the plasma density. V is the velocity. B is 
the magnetic field, p is the plasma pressure, and qo is 
che production rate at a reference altitude h n . H 0 is the 
scale height of the hot neutral oxygen population, ana 
the ratio of speciric heats is 7 = 5/3. 
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In the present first-oraer approximation, the newly 
created photoions are considered as cold (T = 0) and 
motionless (V = 0) particles. A more rigorous approach 
would be to add corresponding source terms to the Eu- 
ler (2) and energy (4) equations ie.g., Gombosi et a/., 
19941. However, it was shown by Biermann et al. [1967] 
that the major effect is that associated with the mass 
continuity equation and that contributions to the mo- 
mentum and energy density are small. Therefore in 
these preliminary calculations source terms in the mo- 
mentum and energy equations are neglected. The ML 
term is given by the right-hand side of the mass continu- 
ity equation (1). Oxygen photoionization is considered 
to be the only mass source with an altitude distribution 
described by the exponential function. This production 
rate is based on PVO observations for solar maximum, 
where the hot neutral oxygen atmosphere had a scale 
height Hq = 400 km in the subsolar region and a refer- 
ence production rate of<j 0 = 3xl0 5 cm" 3 s“ l at altitude 
ho — 400 km \ Belotserkovskn et ah, 1987]. As numerical 
simulations, which were made for various mass addition 
rates, revealed that the boundary layer develops when 
oxygen ion production rate is q 0 = 3 x 10 5 cm* 3 s‘ l 
[McGary, 1993; McGary and Pontius, 1994], we carried 
on our simulations for this value of q 0 . 

Numerical Model 

These numerical calculations were performed with 
HEMIS3D, a three - dimensional ideal MHD code. The 
code utilizes the modified Lax-Wendroff scheme devel- 
oped by Rubin and Burstem [1967], which yields stable 
solutions by adding numerical diffusion to the scheme. 
The dissipative terms (kept as small as possible) are 
included to help damp out short- wavelength ripples 
generated by numerical dispersion and numerically in- 
duced reflections from the simulation boundaries, while 
leaving the longer-waveiength phenomena minimally af- 
fected. The equations are solved in all three spatial di- 
mensions of a spherical (r, 0, d>) coordinate system in 
which the 0 = 0 axis is directed toward the Sun or into 
the solar wind flow. The generally small deviations of 
the solar wind flow from the radial direction (with re- 
spect to the Sun) are neglected so the solar wind flow 
impacting the Venusian ionosphere can be assumed par- 
allel to the Sun-Venus line i0 — 0). The spherical co- 
ordinate system is centered on the planet so r = 
represents the planetary surface. 

Although a grid defined in spherical coordinates in- 
troduces some computational complexity, it has a num- 
ber of advantages for studying flow past planetary ob- 
stacles such as Venus. One obvious advantage is that 
the boundary at the planetary surface occurs at a fixed 
value of the radial coordinate, which avoids the compli- 
cations of interpolating in a Cartesian coordinate sys- 
tem. Another distinct advantage is that as the center 
of Venus is located at the center of the spnerical coor- 
dinate system, high resolution near Venus is obtained 
without the expense oi high resolution evervwnere m 
the simulation. The coae is constructed in such a man- 


ner that it can be applied to a single r. 9 plane wnen 
the phenomenon being studied is cyiindncaily symmet- 
ric about the poles at 0 — 0 and 6 - 180°, as is the case 
when the magnetic field is parallel to the soiar wind 



Figure 1. Plots of the normalized mass density profiles aiomz 0 
- const lines for P = 0 . 1 , p = 0.5, P =10 ana lor me \L 
ML simulations, j. a j the Sun-Venus \d = n line. me 
terminator (0 = 90°) line, ana to the 0 = !$0° line. \M,s 
loading shifts the bow shock farther away from Venus. 
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flow. In this case, the simulation code uses a 73 x 74 grid 
(f, 9 grid dimensions), which corresponds to A9 = 2.5, 
with a radial extent of the simulations of 3 R v . The 
largest Ar = 0.083/Zv. We have checked that the re- 
sults for larger extend of the numerical boundaries were 


very close to those presented in this paper. For the cor- 
responding discussion, see also Cable and Stetnoifson 
[1995]. ' 

The time step limit for methods using explicit tempo- 
ral differencing for advective problems is the weil-known 



Figure 2, 
profiles. 


$ in Figure l out tor the normalized velocity Figure 3. As in Figure i but for 

profiles. 


the normalized pressure 
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Pla ** \ Normalized mass densities {p/ Pao ) for the MHD non loaded (top) and mass loaded 
(bottom) simulations for the case of a parallel (B||V) magnetic field and the plasma 0 = 0.5. 
Here, p* is the density of the solar wind. The effect of mass loading (ML) is to increase the 
mass density m the overall region, except the tail region in which the mass density is reduced. 
Poo denotes the solar wind mass density. 


CFL condition. At < Al/\V\ f where the cell size is A l 
and V is the fastest physical speed in the system. The 
resulting computer code [ Cable and Steinolfson, 1995] 
is second-order accurate in space and time. 

Defining a spherical boundary for the problem is 
easily accomplished. However, solving the equations 
of magnetohydrodynamics in spherical coordinates re- 
quires some modification of the usual Lax-Wendroff 
scheme. In particular, the singular points at 0 ■=. 
0 and 0 — 180° are handled through application of 
L’HospuaTs rule to extrapolate values from grid points 
next to the singular points. Since the value of each vari- 
able at the singular points may vary depending on the 
angle o at which the singular points are approached, a 
flux- weighted average of the values at all angles is used 
tor ttie final updated value. This procedure is described 


in more detail by Cable and Steinolfson (19951. 

We used a simple model for the boundary conditions 
that is physically motivated and performs adequately 
well. Venus is modeled as a hard, highly electrically 
conducting sphere. The inner boundary is taken at the 
planetary surface (1 it*) where the radial components of 
both the magnetic field and velocity are set to zero. All 
other physical variables are computed by extrapolation 
of values along a radial line. At the dayside of the outer 
radial boundary, all quantities are set to the solar wind 
values. The values at the nightside of the outer radial 
boundary are determined by extrapolation along the 
local flow direction. Moreover, we have never encoun- 
tered transient conditions where the flow was coming 
in from the outer radial boundary. A typical compu- 
tation begins with the introduction of the desired solar 
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Plate 2. Normalized plasma flows ( Vfa <») for the MHD non loaded (top) and mass loaded 
(bottom) simulations for the case of a parallel (B||V) magnetic field and the plasma 3 — Q.5. 
Here, a«, is the sound speed of the solar wind. The effect of ML is to decrease the plasma flow. 


wind values in the dayside within the numerical box. 
The initial IMF is taken from the potential solution for 
flow over a sphere so that V B = 0 in the initial state. 
The continued satisfaction of V • B = 0 is assured by 
solving a Poisson equation and updating the magnetic 
field. The numerical solution then continues until the 
interaction achieves an approximate steady state. 

The density and pressure are initially constant through- 
out the simulation region. The initial velocity field is 
the flow of an incompressible fluid over a sphere, di- 
rected along the 6 = 0 line: 


is initialized in a similar way, except that B is skewed 
with respect to the solar wind flow by the desired angle. 
In the present paper, this angle is 0°. This is clearly 
not a frequently occurring case [Phillips et a/., 19861, 
but it is a reasonable first-order approximation which 
will serve as a reference for comparison with results for 
cone angles of 45° and 90°. These processes will be 
realistically modeled in future calculations. 

Results and Discussion 


VJ.(r, d ) 0, t _ Q) _ _ Vq[ 1 — ( cos (6) 

r 

V 9 (r, 8,<p,t = 0) = V' 0 [l + i( — ) 3 ] sin d, { 7 ) 

2 r 

Vi(r, = U) - 0. i 8) 

where Vq is the solar wind speed. The magnetic field 


We present all numerical results for the oxygen ion 
production rate q 0 = 3 x 10 5 cm' 3 s‘ l [McGary and 
Pontius , 1994] and discuss several values of the plasma 
P = 87rpo/i?o» namely 0 = 0.1, 0 = 0.5, 0 = 10, and a 
typical value of 0 = 1.5. The case of 3 = 0.1 (0 = 10) 
corresponds to the strong (weak) magnetic field case. 
The physical solar wind values used are: electron den- 
sity n t - 20 cm' 3 , temperature T = 10 5 °K, sound 
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Plate 3. Normalized plasma pressure (p/(pao a «)) for the MHD non loaded (top) and mass 
loaded (bottom) simulations for the case of a parallel (B||V) magnetic field and the plasma 
0 = 0.5. The effect of ML is to decrease the plasma pressure. 


speed 55 km/s, flow velocity 370 km/s, sonic Mach num- 
ber 6.73, and thermal pressure 6.06 • 10“ l ° dyn/cm 2 . 
Other physical parameters are chosen the same as in 
the work by Cable and Stemoifson [19951. 

As a first step to apply the methodology to the prob- 
lem of ML in the atmosphere of Venus, we consider an 
axisymmetric case with the magnetic field lines parallel 
to the solar wind flow. Consequently, B is parallel to 
the Sun- Venus line and the entire solution is rotation- 
ally symmetric around this line, d/d<t> = 0. We used 
the numerical technique described in the previous sec- 
tion to obtain a solution of the set of (l)-(5). The flow 
is from left to right and all physical quantities are either 
symmetric or antisymmetric about the lower boundary. 

As shown in Plate i, the mass density is increased in 
the case of mass loaded (ML) solutions. This finding 
is consistent with a simple analysis of mass continu- 
ity equation il) and with previously published resuits 


[e.g M Linker et a/., 1989; McGary and Pontius, 1994). 
However, a more detailed close-up reveals that in the 
tail and terminator regions the mass density is reduced 
by the ML effects (Figure 1). Similar reduction is ob- 
served in the case of 0 = 0.5 along the Sun- Venus line, 
at r ~ 1.07 fZ* (Figure la). The density reduction 
does not occur at the terminator for 0 — 0.1, while 
this effect takes place for 0 - 0.5 and 0 = 10 (Fig- 
ure lb). The density depletion is stronger for a higher 
value of 0. Just at the planetary surface, at r < 1.3 
it,, the mass density increase is observed (Figure lb). 
In the magnetotail, the plasma density reduction occurs 
at the planetary surface for all values of the plasma J. 
Again, the rate of density depletion is proportional to 
0. The case of 0 = 0.1 differs from the other cases just 
as at the planetary surface p is increased, but farther 
on there is a region of mass depletion which is followed 
by a region of mass enhancement (Figure lc). To our 
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Piate 4. Normalized plasma temperature (T/T.) for the MHD non loaded (top) and mass 

/? _ ne simulatlons for case of a parallel (B||V) magnetic field and the plasma 
p — 0.5, The effect of ML is to decrease the plasma temperature* 


knowledge, this is the first report of mass depletion due 
to the ML. The density depletion is a consequence of 
the nonlinear term V ■ (pV), which is smaller due to 
the plasma flow reduction by the ML effects. Indeed, 
the largest density (velocity) increase (decrease) occurs 
in the magentotaii at r = 3 for 3 = 0.1. The nor- 
malized density jumps from p NL - 1.25 to p ML - 1.5 
(Figure ic), while the magnitude of plasma flow is re- 
duced from V NL cr 4.9 to V ML - 4 (Figure 3c). Con- 
sequently, pv L V NL = 6. 125 > pmlVml = 6. Therefore 
the effect associated with this nonlinear term overcomes 
the density increasing effect exerted by the exponential 
term go exp (-(h — ho)/ Ho]. The density reduction in 
the tail may also be due to the changed shape of the 
bow shock. A different bow shock shape may result in 
less mass being carried into the tail, regardless of other 
factors. 

It follows trom nuler equation ( 2) that we can expect 


the reduction of a velocity magnitude because conserva- 
tion of momentum and energy requires that the plasma 
flow is decelerated as a consequence of the mass density 
reduction. Plate 2 presents results both for the ML and 
no loading (NL) simulations. By examining the two so- 
lutions, one can see a fundamental difference between 
them. The effect of ML is to decrease the flow. A de- 
creasing rate is highest in the tail region where slow 
flows in the case of NL are much more reduced by ML. 
The decreasing rate is larger for a smaller value of 3 at 
the terminator and in the tail region. The plasma starts 
decelerating because of ML well ahead of the shock. In 
the subsolar region, this rate is largest for 3 - 0.5. A 
large flow deceleration occurs also at the terminator, 
close to the planetary surface where the plasma flow 
drops from V SL ~ 5.7 to V ML ~ 4.2 (Figure 2b). In- 
spection of Plate 2 reveals the global structure of the 
shock and the stagnation region behind it. The flow 
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decelerates at the shock as approaching solar wind ab- 
sorbs an increasing number of ions. At the planetary 
surface, r = there is no perpendicular flow as a con- 
sequence of the boundary conditions imposed. Velocity 
shears in the radial direction are present at the bow 
shock region both for the ML and NL flows, which is 



r 


Figure 4. As in Figure i but for the normalized temperature 
profiles. 


consistent with PVO observations iMihalov et ai % 19821 
and the g&sdyn&mic results by McGary [1993], 

From (4) we can derive an equation which governs 
the evolution of the NL plasma pressure p 

dp 

— + V.(pV) = ( I- 7 )pV.V. (9) 




Figure 5. As in Figure 1 but for the normalized magnitude of the 
magnetic field profiles. 
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As the plasma flow velocity V is decreased by the ML 
effects we can expect a pressure reduction. The simula- 
tion results shown in Plate 3 are qualitatively consistent 
with our expectations. However, Plate 3 ^.hows that the 
plasma pressure is increased near the terminator region, 
in the region 9 < 90°. This increase is caused by a 
nonlinear interaction which is described by the V (pV) 
and pV* V terms. Our results differ from the conclusion 
made by McGary and Pontius [1994], who claimed that 
the dynamics of the plasma flow are not significantly 
affected by ML. A reason for the difference between the 
present results and those of McGary and Pontius has 
to do with the different models applied; McGary and 
Pontius performed gasdynamic (with B = 0) simula- 
tions for the cylinder flow, whereas our model concerns 
MHD and a plasma sphere. 

The above scenario may differ for different values of 
/3. It occurs in the case of a sufficiently small value of 0 
that the magnetic field is so strong that the fast shock 
becomes an intermediate one [ Stetnolfson and Hund- 
hausen, 19901. For example, the ML increases the gas 
pressure along the whole Sun- Venus line in the case of 
d = 0.1 and 0 = 10 (Figure 3a) but reduces it in the 
case of 3 = 0.5, 0 = l, and 0 = 1.5 at the planetary 
surface. In the bow shock region, an enhancement of 
p due to the ML is observed for all values of 0. The 
strongest enhancement seems to occur for 0 = 0.5. At 
the terminator (Figure 3b), the gas pressure is increased 
by ML. The rate of decrease is highest at the planetary 
surface for 0 — 0.1, while at larger distances the influ- 
ence of ML on p is negligible. The pressure reduction at 
the planetary surface is caused by the right-hand side 
term of (9) which becomes small at the surface because 
of the flow reduction. The bow shock is closer to the 
planetary surface for a higher value of 0, in accordance 
with our expectations. At the magnetotail the pressure 
is reduced (Figure 3c). The reduction rate is highest at 
the planetary surface for the smallest value of 3 and is 
a monotonic function of 3. 

The next point of comparison between the ML and 
NL results is the temperature distribution around Venus. 
As for the ideal gas the plasma temperature T is pro- 
portional to p/p so both p and p play a major role in 
determining its behavior. We already learned that for 
0 = 0.5 p (p) is increased (decreased) by the ML effects 
everywhere except in the terminator and tail (nose) re- 
gions. Therefore, we may expect the general behavior of 
T would be to decrease its magnitude by the ML effects. 
Indeed. Plate 4 indicates that the ML significantly cools 
the flow relative to the NL case. The cooling occurs as 
zero-temperature oxygen ions are introduced into the 
flow and thermalized [McGary and Pontius . 1994]. The 
behavior of T is interesting to discuss in more detail as 
a function of the radial direction r at different locations 
of 9. In the case of 9 = 0. for d - 0.1 and 3 = 10 
the temperature is decreased by ML both at the plane- 
tary surlace and at the sunward side of the bow shock 
(Figure 4a), while in the shock the effect of ML is neg- 
ligible. Tne curves which correspond to the ML possess 
maxima in the magnetosheatn. wmie no maxima exist 
tor the NL case. A similar maximum occurs also in the 
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Figure 6. Terminator positions for different values of the piasma 
P for the non loaded and mass loaded simulations. The bow 
shock at the terminator is shifter farther away farther away from 
the surface of Venus. 

case of 0 = 0.5 for the ML curve. Now, however, ML 
cools the plasma at the planetary surface but warms it 
in the bow shock region. The temperature increase at 
the bow shock was observed recently by Gombosi et ai 
[1994] in their simulations of the interaction of an ex- 
panding cometary atmosphere with the solar wind. At 
the terminator, the plasma is warmed by ML (Figure 
4b). The warming is more significant at the planetary 
surface and is higher for a higher value of 0. The cooling 
occurs at the magnetotail (Figure 4c). The cooling rate 
is higher at the surface of Venus and for a smaller value 
of 0. This is a consequence of the pressure decrease 
at the planetary surface (Figure 3c). The mass den- 
sity is increased for 0 = 0.1 and decreased for 3 = 0.5 
and 0 = 10, but the decreasing rate is smaller than 
the pressure decreasing rate (Figure lc). Consequently. 

T ^ p/p is much decreased at the surface of Venus. The 
cooling rate is much more significant than the warming 
effect. The bow shock position is closest to the plane- 
tary surface for the smallest presented values of 3. The 
bow shock distance from Venus grows with 0. 

Plate 5 shows the magnitude of the magnetic field. 

In the subsolar region there is no jump in the magni- 
tude of the magnetic field across the bow shock. At the 
flanks, however, the magnetic field increases through 
the bow shock. The magnetic tail is a region of reduced 
magnetic field. It can be seen that the ML exerts some 
influence on the magnetic field strength. Small differ- 
ences occur at the sunward side of the shock where the 
magnetic field is decreased by the ML effects (Figure 
5a). This effect is small for 0 — 0.1 and grows with 
0. The decrease in the magnetic field strength nearest 
the planet along the stagnation streamline, is the re- 
sult of the boundary conditions (B r (r - R v ,6 y t) ~ 0) 
imposed together with the fact that the parallel shock 
does not change a value of the magnetic field across 
the shock. As the planet is conducting there can be 
no radial magnetic field component at the surface oi 
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Venus. But near the subsolar line, the solar wind mag- 
netic field is purely radial. Consequently, the magnetic 
field strength has to drop to zero in this region. This 
behavior counter to the common picture of B piling up 
in the stagnation region in the case of B 1 V [Ca6/e 
and Stemolfson , 1995; Murawski and Steinoifson , 1995] 
and would only occur for the magnetic field nearly par- 
allel to V. At the terminator, the ML enhances the 
magnitude of B at the bow shock. Therefore at the 
planetary surface B is reduced. A decreasing rate is 
hardly dependent on 0. The largest differences occur in 
the magnetotail where the magnetic field is enhanced 
by ML. The strongest enhancement occurs for /3 = 0.1. 
This behavior contradicts our expectations, which fol- 
low a simple analysis of (3) wherein as a consequence of 
V decrease we expect the magnetic field reduction by 
the ML effects. However, V couples nonlinearly with 
the magnetic field. The coupling occurs through the 
term V x (V x B), which determines the overall behav- 
ior of the magnetic field rather than the plasma flow V 
only. 

Inspection of Plates i-5 reveals that a bow shock is 
formed upstream of Venus. The principal effects are 
that the bow shock is displaced farther from Venus in 
the case of the ML plasma [Belotserkovskii et a/., 1987; 
Breus et a/., 1989; McGary and Pontius , 1994]. For 
precise investigation of the bow shock structure, Fig- 
ures 1-5 exhibit the plasma mass density, velocity, gas 
pressure, temperature, and magnetic field distribution 
along the Sun-Venus [6 = 0°), terminator (6 = 90°), 
and 0 = 180° line. Inspection of these figures reveals 
that the pressure at the planetary surface is larger for a 
smaller value o( 0. The ML effects also shifts the shock 
at the terminator farther away from the planetary sur- 
face (Figure 6). This shift is largest for intermediate 
values of the plasma 0. These results show the bow 
shock at the terminator is in about the position of the 
PVO observations [Tatraliyay et ah. 1984: Russell et 
al.< 19901. Average subsolar ana terminator distances 
of the bow shock are about 1.15 R v and 2 R, (from 
the center of the planet), respectively. To match the 
observed locations (at the nose) of the bow shock, the 
MHD model requires that the mass loading be greatly 
enhanced. 

Russell and Zhang [1992] in their analysis of PVO 
data found that during periods when the solar wind 
magnetosonic Mach number is near unity and the plasma 
beta is low, the bow shock may travel to distances 
greater than 10 R v away from the planet. On the ba- 
sis of previous work on the formation of MHD shocks 
in coronal mass ejections [Steinoifson and Hundhausen , 
1990], it was felt that the solar wind conditions had be* 
Lome {during the far bow shock excursions) such that 
the usual fast MHD shock is no longer a stable solu- 
tion. The lower than normal plasma beta and mag- 
netosonic Mach number are in a parameter regime for 
which the usual fast-mode bow shock close to the planet 
may not provide the necessary compression and deflec- 
tion ol the solar wind. Lsing MHD simulations, it was 
snown ! Steinoifson and Cable, 1993: that, for these con- 
ditions. tne usual fast shock is repiacea by a Dow shock 



Figure 7, Time history of the radial velocity component ( l V i 
at the spatial position very close to the bow shock for the non 
loaded simulations. The upper ( lower) curve corresponds to d 
10 (P =0.1). Not the stabilizing effect of the stronger matinc::. 
field. 

configuration containing an intermediate shock near the 
Sun-Venus line and a fast shock at large distances from 
the Sun-Venus line. The resulting configuration propa- 
gates upstream away from the planet at a low speed and 
appears to be approaching a new equilibrium stand-off 
location at a large distance from the planet. 

Plates 1-5 exhibit waves which propagate along the 
bow shock boundaries. We have checked that these 
waves are not numerical artifacts as they do not dis- 
appear for a finer numerical grid and are persistent in 
time (not shown). These waves are generated by the 
solar wind impact, which is exerted by the solar wind 
on the plasma surrounding the planet. The first impact 
shifts the plasma towards the surface of Venus, while 
later in time the magnetic field recovers and the plasma 
is displaced farther away from the planet. Moreover, 
the bow shock first builds closer to the planet and then 
its altitudes at the equator (6 - 0°) and at the termi- 
nator (0 = 90°) grow in time. As a consequence, the 
nose of the bow shock oscillates in time. The Alfven 
and slow oscillations, which are guided by the magnetic 
field lines (e.g., Murawski et a/., 1995], are fed by the 
Kelvin-Helmholtz (KH) instabilities [e.g., Rankin et a/.. 
1993], which are more apparent in the case of a weak 
magnetic field case. Figure 7 presents the radial com- 
ponent of the plasma flow as a function of time and for 
3 = 0.1 and 0 = 10. This component is calculated at 
a spatial position which is very close to the bow shock. 

As expected, the stronger magnetic field makes the flow 
less unstable as the amplitude of oscillations is smaller 
in the case of 0 = 0.1 than in the case of 3 ~ 10. The 
KH instabilities criterion [e.g., Rankin et a/., 19931 
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( Bo * k r 


'.Pi + P2T 


4*(Pi t- f j 2 ) 


( 10 ) 



2558 


MURAWSKI AND STEINOLFSON: SIMULATIONS OF MASS LOADING 




Plate 5. Normalized magnitude of the magnetic field (B/y/p opooajj for the MHD non loaded 
(top) and mass loaded (bottom) simulations for the case of a parallel (B||V) magnetic field and 
the plasma 0 - 0.5. The effect of ML is to increase the magnitude of 5 in the terminator region. 


is satisfied easier for a less magnetized flow. In this for- 
mula indices 1 and 2 correspond to variables on different 
sides of the bow shock, while denotes the growth rate 
of the KH instabilities. It is easy to find out that the 
flow is unstable (w? > 0) for a non-magnetic plasma 
(Bj = B 2 = 0). On the other hand, strong magnetic 
fields can make the flow stable (u;? < 0). Intermediate 
values ot the magnetic fields stabilize the flow and in 
particular the flow can be marginally stable (u? = 0). 
These instabilities occur in the region where the flow 
has a large gradient such as at the bow shock. Assum- 
ing that the wave propagation k is parallel to the flow 
V and magnetic field B, it follows from (10) that 


(Pi +P7 = PxPi(V 2 - V. r - ■ /n -rp 2 )(Bl + Bj), 


ill) 

where now the plasma parameters are normalized by 


the units of the solar wind. V A is the Alfven speed of 
the solar wind. At the surface of Venus , Vi — Bi = 
0. Upstream of the bow shock, the mass density p l 
and the plasma flow V\ are hardly dependent on the 
plasma 0. Therefore we can put p x ~ 0.95 and V\ ^ 6.7 
in normalized units. For 0 = 0.1 (0 = 10) we have 
P2 ^ 3.9 (p 2 - 3.7) and B x ~ 0.93 (B x cr 0.91). It is 
easy to check that the right-hand side of ( 1 1 ) is positive, 
indicating that the flow is KH unstable for both 0 - 0. 1 
and 0 = 10. 

Summary 

The interaction of the solar wind with Venus is stud- 
ied using numerical solutions of the time-dependent 
ideal two - dimensional MHD equations. A modified 
Lajc-Wendroff scheme [Cable and Steinoifson, 19951 is 
used to solve the equations in a spherical coordinate 
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system. For these computations the detailed chemistry 
of the ionosphere is neglected, and the planet is treated 
as a conducting sphere. Numerical computations have 
been perfomed to study the effect of ML and the IMF 
strength on the magnetic barrier and the general con- 
figuration of the magnetic tail. The IMF orientation 
included in this study is parallel to the solar wind flow 
direction. 

The mam results are the following: The ML increases 
the mass density in the overall region, except the termi- 
nator and tail where the mass density depletion occurs 
as a consequence of the nonlinearity action. The mass 
density depletion is larger for a higher i3. As a conse- 
qence of this reduction, the plasma flow, pressure, and 
temperature are reduced. This scenario depends, how- 
ever, on a value of the plasma fi. It was found that for 
0 — 0.1 both the plasma pressure and the temperature 
are increased at the dayside by the ML. 
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Abstract. The large-scale (global} interaction of the 
atmosphere of Venus with the surrounding medium is 
numerically simulated in the framework of 3-dimen- 
sionai magnetohydrodynamics. The simulations are 
performed for the case of the interplanetary magnetic 
field perpendicular to the solar wind flow. The model 
reproduces several features of the interaction predicted 
by earlier theories and observations. The soiar wind 
interaction with Venus is charactenzed by the mass 
loading of the solar wind with heavy oxygen ions which 
are produced by the photoionization of neutrals in the 
extensive ionosphere. This mass loading slows down 
the solar wind and ultimately leads to the formation of 
a magnetic barrier and a magnetotail. The location of 
the bow shock is barely affected by different values of 
the magnetic field strength. 


I. Introduction 

The soiar wind is an almost coilisioniess plasma consisting 
mainly ol protons and electrons which flow outward from 
the Sun, eventually exerting an impact on planets. 
Although the solar wind carries with it the interplanetary 
magnetic field (IMF), Venus has virtually no intrinsic 
magnetic field (e.g. Luhmann. 1986) Direct observations 
< by the Manner. I enera. and Pioneer lenn.s missions) at 
Venus have demonstrated that the planet's atmosphere 
consists of a magnetosheath bounded b\ a bow shock on 
the outward side and ionopause on the lower side. The 
bow shock is a thin region across which the plasma par- 
ameters change drastically. The bow shock exists because 
the solar wind is supersonic and super-Allvenic and must 
become subsonic in order to flow around the obstacle. 
The ionospheric plasma which tries to exclude the solar 
wind plasma with its associated IMF is a verv good elec- 
trical conductor Hence, the tonospnere forms a nearly 
impenetrable obstacle which can serve to divert the solar 
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wind around it. The shocked solar wind slows down a> o 
approaches the planet, and as a consequence, the muanem 
held strength builds up. resulting in the formation o; 
magnetic barrier located outside the dayside ionopause 
The ionopause is usually coincident with the bottom m 
the magnetic barrier. It is one of the major aims of space 
plasma physics to measure the plasma properties through- 
out the whole atmosphere and to develop computer mod- 
els for its representation. Evaluation of a computer model 
must depend ultimately on how- well the results match 
observations, but valuable insight can also be gained b\ 
comparison of the results that were obtained by different 
numerical techniques and by analytical estimation ii av ail- 
able. 

Application of magnetohydrodynamics to the soiar 
wind interaction with Venus was initiated by Sprciter et 
aL (1970). At that time the gasdvnamic equations were 
solved for the density, velocity, and pressure, while the 
magnetic field was evaluated in a subsequent step b> sow - 
mg the induction equation using values for the densitv 
provided by the gasdvnamic solution. This approach con- 
cerns the scenario in which the magnetic field lines arc 
convected with the local flow of the plasma. 

Mass loading (ML) effects have been implemented b\ 
Belotserkovskiie/ al. (1987), Cloutier etai. (1987), Phillips 
et al. (1987), Stahara et al. (1987), Breus et al. ( 1989i. ana 
Sauer et al. (1994). A comparison of the results of Stahara 
et al. with those of Belotserkovsku et al. reveals several 
major differences. The results of Belotserkovsku a al. 
indicate that the ML has a large effect on the location ot 
the bow shock, moving it to a position slightly bevonu 
that observed by PVO, while Stahara et al. claim that me 
location of the bow shock is barely affected by the ML oi 
the flow. 

It has been shown by Murawski and Sleinoifson i mviM 
that in the case where the IMF is parallel to the solar w mo 
(low. the ML increases the mass densitv m the overall 
region, except the terminator and tail where me mass 
density depletion occurs because of the nomineuntv 
action. As a consequence of this reaucuon the niasnnt 
fiow. pressure, and temperature are reaucea. Tins scenario 
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ueoencis. however, on a value or the masma /. It was 

I'm ?, 0 1 " ll 101 ~ '' ' lde h'asma pressure and the 

vmpe rat ure are increased at the dayside by the ML. 

^ f IS 1 ie P UI P° se this paper to test the conclusions of 
a ara (/(/. and Belotserkovskn et ai.. and extend the 
analvsts oi Muraw sk. and Ste.noitson into ihe case of the 
1 ,h Perpendicular to the solar wind flow 

The paper is organized as follows. Mathematical and 
uumciical models ol the dynamics oi a venusian plasma 
' r«; S1 ?'\ n 111 '^ CCII ° n • s,umer| cail> obtained results are 
P vnh C ? m a UlCr SeCtl0n ol this P J P er - The paper closes 

,th Dnet sumrn ary of presented results and general 
remarks. 


2. Simulation model 

The interaction between a solar wind and the atmosphere 
enas ls as sumed to be governed bv ideal MHD equa- 
tions which represent coupling of the fluid dynamic equa- 
lonswith Maxwell s equations of electrodynamics by nee- 
-sting electrostatic force, displacement current heat 
■'-ansier. and damping etfecis We solve the lollowimi form 

problem° mPrCSSlble MH ° Cc ! UUUons ^ an imtiafvalue 


K ».»•, ana R. S S«moif S on : Num t „cai „„„ y<|> 


cp 

-rj t-VMpV) = \l 


(I) 


U> V) I 

•— -V-[ (/ A)V] =-V/>+-iVxB|xB 1 2 ) 


Ct 




l> V 


+ -[Bx|VxB)J 


cB 

ct 


— V x ( V x B ) 


0 (3) 


(4) 


tvl re ,u.'n the m aS densily ’ , V the velocity. B the divergence 
0) magnetic field, p the gas pressure, and the 
ratio ot specific heats is y = 5 '3. 

We assume that the only important source of ion pickup 
is due to the photoionization of the hot oxygen corona 
Ihe newly created ions move so slowly thaT essentially 
their flow does not contribute to momentum and energy 
equations. Moreover, we consider the case that the ion 
production rate is proportional to neutral atmosphere 
oxygen density In the present first-order approximation, 
the newly created photoions are considered as cold 
■ M and motionless |V = 0) particles. A more rig- 
orous uay would be to add corresponding source terms 

° E | U ooav , 'L and energ> <3) e P uatl °ns (e g. Gombosi et 
;' 0 ' 99 y However, u was shown by Biermann et at. 

( I %7) that the major effect is that associated with the mass 
continuity equation and contributions to the momentum, 
magnetic field, and energy density are small. Therefore. 
111 r ^ ese preliminary calculations source terms in the 
momentum, magnetic field, and energy equations are neg- 


; ecied. The mass ioading term M can tnus oc represent 
as (McGary and Pontius. 1994) 

•V/= M 0 exp ;,//„] , 

where M is proportional to the product ol the numn 
density ol the neutral hot oxygen molecules ai the plane! 
nose and their mass, divided by the characteristic mne i, 
photoiomzation. The magnitude of the mtensiiv oi >, 
Pickup at Venus appears to vary significantly over ihe 
vear solar cycle (Russell et al. % 1990). 

The ion production rate is based on Pioneer ( ,/,,■ 
observations for solar maximum, where the hot neutr 
oxygen atmosphere had a scale height //,, = 400 km i 
the subsolar region and a reference production raic < 
•V/« 3 x10 cm s at altitude /i„ = 400 km (Iklo 

serkovskn at., 1987). As numerical simulations mac 
for vanous mass addition rates revealed that the bounaar 
layei^ develops when the oxygen ion production raic 

m.° Pools 0 Cm 5 (McGa ry, 1993: McGarv and Pui 
bus 1994), we earned on our simulations for this vah, 

OI Aj q . 

The MHD equations are solved numerically using in 
modified Lax-Wendroff algorithm (Rubin and Buimc.i 
, 7 >; ! n our 3 -dimensional s.mulations the equations ai 

;? V ® d 'I 3 sphencal (r - 9 ' <£) coordinate system in whic: 
the 6 -0 axis is directed toward the Sun and 0 = 0 cor 
responds to the equatorial plane. This coordinate svsten 
is centered on the planet which is modeled as a hard 
highly electrically conducting sphere. Consequently, ih 
plasma velocity and magnetic field arc required to b. 
parallel to the surface of the planet, and zero-order rauia 
extrapolation is used to determine the other quantities oi 
the planet surface. The values at the mghtside of the outci 
radia! boundary are determined by extrapolation alone 
the local flow direction. To carry out numerical siirui- 
lations, solar wind conditions are maintained on the dav- 
side of the outer radial boundary. The simulation code 
used a 73 x 74 x41 grid (r. d. 0 grid dimensions). Thb 
corresponds to the angular grid spacines. AU = ■ - ,, 1U 
A0 = 4.75, with a radial extent of the simulations of 3 ft, 
The largest radial spacing is Ar = 0.083 R v . The resulting 
computer code (Cable and Steinolfson. 1995) is second- 
order accurate in space and time. 

The denisty and pressure are initially constant through- 
out the simulation region. The initial velocity held is the 
flow of an incompressible fluid over a sphere directed 
along the 0 = 0 line : 


K(r, 0, 0, t = 0) = V, 


K)(r, 0, 0, / = 0) = W„ 


/?, 


1 //?v 


- \ r 


COS U 


sin u 




K,(r,0,0,f = 0) = o. 

The magnetic field is initialized in a similar wav. except 
that B is skewed with respect to the solar wind flow bv tn ■■ 
desired angle. The initial IMF is taken from the potential 
volution for flow over a sphere so that V - B = ii m me 
initial state. The continued satisfaction oi V-B = .« 
assured by solving a Poisson equation ana updating me 
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Hg. I. Overall behavior of t he solar wind- Venus interaction for the case of the perpendicular magnetic 
field. The initial locations of the blue fluid elements were equally spaced along a line located well 
ahead of the shock front, perpenaicuiar to the solar wind flow. Their positions are shown in the heure 
at tune increments ol 5 s. Note that the fluid elements are slowed as they entered the shock front, then 
accelerated again when they passed over the planet surface. The magnetic field lines, in ycilow. piied 
up at the planet surface 



Hi-. 2. Reconnected lieid hno in the mil region in the case ot the magnetic held perpendicular to the 
oiar wind how ,ne \L (b) the ML case. Note that as a consequence of the plasma riow 
reduction, the reconnected linos are pulled away trom tne planet dv the ML effects 
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^oiur wind interacuon uun Venus 



My, 3. A summary of ihe lie Id line behaviors m the NL (a) and ML (b) eases. The view is from i he 
tail looking at the mghiside of the planet 
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\Ls m It i ,l SS !h 1 ! ilri i ll ' IOm °! thC SOk,r W ' nd : (U) the e chptic and noon-m.dnight planes for the 

nlan^for the \ I uf'”" 0 noon-m,dn,ght planes tor the ML s.mutat.ons : (c) the terminator 
niane lor me .\L and ML simulations 
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oagnetic nciu. The numerical 'imiiion men continues 
until the interaction acmes es an approximate steady state. 
\ typical computation begins sum the introduction of 
me desired solar wind \aiues in the day side within the 
numerical box. 


3. Numerical results 

We present aii numerical results tor the oxygen ion pro- 
duction rate A/,, = 3x 1 O' cm > ( McGary and Pontius. 

1994) ana perform computations with a cone angle of 90 
»the perpendicular case), while the average solar wind 
cone angle is reported as 42:2 (Phillips a ai.. 19X6). The 
computations wuh the cone angic of 0 (the parallel case) 
have been reported previously by Murawski and Ste- 
moitson ( 1995). The physical solar wind values used are; 
electron density n, - 20 cm . temperature T = UP K. 
>ound speed 55 km s', flow velocity 370 km s' ! . sonic 
Mach number 6.73. and thermal pressure 6.06 x 10 “ 10 dyn 
-m ' Other pnysical parameters arc chosen the same as 
a Cable ana Steinoifson i i49u 

h mure l shows the global connguration of the magnetic 
held lines and fluid elements traces from the final con- 
figuration of the numerical simulation. The solar wind 
hows m Irom the left-hand side towards the planet. Yellow 
lines indicate magnetic held lines which pile up at the 
shock tront. and then slip over the planet, forming mag- 
netotail. The blue dots represent the instantaneous pos- 
itions oi nine fluid elements moving from left to right in 
the plate. Successive positions for each fluid element are 
separated from one another by a time increment of 5 s. 
The huid elements that enter the bow shock center are 
significantly slowed in their courses, then accelerated 
again as they pass over the surface of the planet. In the 
davside part of the planet, how velocities become small, 
and a stagnation point is formed ahead of the planet. This 
scenario is hardly affected bv the M L effects and is similar 
ror different vaiues oi the soiar wind magnetic held 
strength. 

Figure 2 shows the reconnected magnetic held lines in 
the magnetotail of Venus. These lines are stretched out 
into the tail in a roughly conical volume with a base radius 
which depends on the strength of' ML. As the ML effects 
reduce a plasma flow in the tail region, the reconnected 
region is shifted farther out from the planet surface for 
non-zero ML effects. The ionosphere of Venus captures 
held lines which drape over the planet to form an anti- 
parallel magnetic held connguration in the magnetotail. 
Converging hows in the magnetic lobes squeeze the anti- 
parallel magnetic lines to form a thin plasma sheet. 
Eventually, magnetic reconnection occurs in this region. 
The resistivity producing the reconnection originates in 
numerical diffusion which inherently builds up in a 
numerical scheme. This scenario resembles reconnection 
events m the magnetotaii of cornets (Ogino cl aL 1986). 
Our findings of the magnetic reconnection in the tail 
region are consistent with McComas ct ai .' s ( 1986) resuit. 
Mat a- 12 R v downstream, tne crossflow no id is reversed 
>mv .mom : oi me time, foe present results suggest 

oso that a mgner vaiuc oi me ML rate is needed to smft 


he reconnection region farther out from tne surface i 
Venus. 

Figure 3 summarizes the various held lines we nave 
observed in the case of no loading (NL) ana ML .emu- 
lations. The held line passing around the opposite side oi 
Venus has flowed downstream with the impacting soiar 
wind and has become caught on the davside of the pianet. 
The uppermost held line slips over the planet s surface 
and bends into a U-shaped form, with the U-openinu 
in the downstream direction. Spiraling held lines winch 
exhibit reconnection can be seen toward the bottom. The 
ML affects mostly the magnetic held lines at the mag- 
netotail (Fig. 3b) where two lines already underwent me 
reconnection, and the third one, in a U-shaped form. i> 
close to a reconnection state. A plasma how in the euua- 
tonal plane brings the two magnetic held lines into the l - 
shaped form. These lines may undergo reconnection if the 
how in the d> = 0~ plane is small enough not to reduce the 
curvature of magnetic held lines. As ML decreases, the 
how in the tail region of the reconnection occurs farther 
out from the pianet surface. 

Figure 4 shows the influence of the ML effects on me 
radial (or B.) ana tangential (or B v ) components oi me 
magnetic held. These components are measured amne 
the circular orbit of r = \AR V in the equatorial [0 = u i 
plane. The values of the radial component are increased 
by the ML effects and are very similar in magnitude to 
those recorded by PVO (Luhmann et ai .. 1991 ; Dubinin 
et ai. . 1991). Although the ML enhances the tangentiai 
component of the magnetic held, its value is still too small 
at the current sheet which is located at 0 = 180 Luhmann 
e/ ai. (1991) show that this component reaches a value of 
^ 10 nT at the current sheet. This disagreement indicates 
that the rate of the ML effects is too small or other pro- 
cesses apart from MHD phenomena are at work main- 
taining the cross-flow held strength in this region. 

Figure 5 exhibits the spatial distributions of the gas 



theta (dea) 


Fig. 4. Radial or r and 0 or .v components of the magnetic ::eui 
in a circular orbit in the ecliptic plane of r = I 4/? v Note mat 
the ML effects increase the magnitude of B in the overall reeion 
of 0. The B x component is increased by the ML effects eve:", 
where except 9^0 (sunward) direction, Note mat me 
numerical simulations 0 <: 9 $ 180. To draw cas net; re 
ne symmetry 'ounaary conditions a ere .mmiea 
SO' < ** < 360 
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Fig. 6. The plasma ip.T.p) profiles in the case of the perpendicular magnetic field for the NL (broken 
lines) and ML (solid lines) simulations. Profiles are made along: (a) the Sun-Venus line: (b) 6 = 90 . 
4> = 90 line ; tc) 0 = 90.0 = 0 line ; (d) 0 = 1 80° line. The grid points were denoted by stars in the 
pronle ol the ML piasma pressure p (a). The mass density p is denoted here by d 


pressure from the perpendicular simulations. As stated 
above, the outer computational boundary was 3/? v . but 
here we show r only the region contained within 1.75/? v in 
order to bring more details. Figure 4a i4b) shows the 
ecliptic and noon-midnight planes for the case of NL 
(ML) simulations. Formation of the bow shock can be 
seen in front ol Venus. The distance between the shock 
and Venus is larger in the case ol ML. This is in agreement 
with observations by PVO. When the solar activity was 
near its maximum these observations indicated that the 
bow shock receded significantly from the surface of Venus 
(e g. Slavin el a/., 1983). As the time progressed toward 
solar minimum, the shock approach ed the planet's 
surface, and the process eventually repeated following a 
trend similar to that ol the previous solar cycie (Tatrailvav 
t ft a 1984). 

Note that the pressure maximum occurs downstream 
ol the shock, then drops to a lower value at the planet's 
surface (Fig. 6a). See also Tanaka (1993). and Cable and 
Stcinoifson (1995). In the parallel case the pressure 
maximum occurs at the planet surface tCable and Ste- 
inolfson. 1995: Murawski and Stemotfson. 1995). The 
effect ol ML is to increase i decrease i the pressure 
upstream (downstream) of the NL bow snock. This 
oehavior is indepenaent of the magnetic neia strenetn. In 


the parallel case a similar gas pressure reduction due to 
the ML occurs for fi = 0.5. For 0 = 0.1 and 10. the gas 
pressure is enhanced by the ML (Murawski ana Stei- 
nolfson, 1 995). Similar behavior occurs for the plasma 
temperature which is higher at the subsolar side of the N 1 . 
bow shock. Such plasma warming is observed also in the 
case of the parallel simulations for fi = 0.5. However, the 
cases of /? = 0.1 and 10 differ from this behavior (Mur- 
awski and Steinolfson, 1995). The mass densitv ^ 
increased at the nose region ( Fig. 6a), as in the case of the 
parallel simulations (Murawski and Steinolfson. 1995), 

Figure 5c shows the pressure distribution in the ter- 
minator plane (0 = 90°). Again, the pressure is increased 
(together with T and p increase), and the bow' shock 
distance from the planet's surface is increased when the 
ML is applied (Fig. 6b). The increase of p occurs mamlv 
in the magnetosheath. Note also that the bow shock is 
thicker in the case of ML. There are two small depletions 
of p which occur for the ML solutions at both the equa- 
torial and day-night planes ( Figs 5 and 6b. c). 

It appears with some effort on Figs 5c and 6b. c mat 
the bow shock is flattened at the terminator plane. The 
flattening occurs in the equatorial plane as the shock dis- 
tance from the pianet surface is larger in the & = u piane 
than in the o = 90 ' plane. Although this is in agreement 
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Fig. 7. Ah tor Fig. 2 but tor 8 and V Profiles are along the Sun- 
Venus line i a ) and along the 0 = 180 line (b) 

with the results by Russell et al. ( 1 990), it is rather unlikely 
that the degree of flattening of the Venus bow shock could 
be accounted for in this wav to match the observations. 
These results are also in agreement with the conclusion 
drawn by Russeil a at. ( 1988) who showed that the cross 
section ot the bow shock in the terminator plane was 
almost perfectly circular when the IMF was aligned with 
the soiar wand flow but departed from circularity by about 
10% when the IMF was perpendicular to the solar wind 
direction. For the asymmetric case the smaller diameter 
was found to be aligned with the IMF 

In the magnetotail region the mass density is enhanced 
at the planet surface due to the ML. At farther distances, 
the ML reduces p (Fig. bd). This behavior differs from 
the parallel case lor which p is reduced at the pianet and 
enhanced in the tail, at larger distances. In the parallel 
case o is enhanced for (1 = 0. 1 but. depletion regions exist 
lor is = ofand 10 I Murawski and Steinoifson. 1995). The 
gas pressure and temperature increase due to tne M L ( Fig. 6d). 

The magnetic held strength rises monoionicaily to a 
peak at the planet's surface indicating me formation of a 
magnetic barrier (Cable and Steinoifson, i995>. A detailed 
investigation of the magnetic held along the Sun-Venus 
; me shows that B is enhanced (reduced) by ML upstream 
■■.:own>ireami oi the non-iouded bow shock iF;g. 7a). At 
• •e terminator m ) B is increased m me ML effects 


while in the parallel case. B is ennanccu at distances farther 
(closer) to the planet’s surface (Murawski ana Stei- 
noifson. 1995). The departure of the bow shock from the 
planet's surface as a consequence of the ML applied is 
also discernible on the spatial distribution of the magnetic 
field. In the magnetotail. ML enhances B at the pianei 
surface (r < 2.5 7%) and reduces it at longer distances ( Fig. 
7b). Note also two minima which correspond to recon- 
nection regions. ML shifts the reconnection rcaion from 
r - 2.357% to 2.87%. 


4. Summary 

The interaction of the soiar wind with Venus is studied 
using numerical solutions of the time-dependent ideal 3- 
dimensional MHD equations. A modified Lax-Wendrotf 
scheme (Cable and Steinoifson. 1995) is used to solve 
the equations in a spherical coordinate system. For these 
computations the detailed chemistry of the ionosphere is 
neglected, and the planet is treated as a conducting sphere 
Numerical computations have been performed to stuuv 
the effect of ML, the IMF orientation and strength on the 
magnetic barrier, draping of the magnetic field over the 
planet, the slippage of the field around Venus, and the 
general configuration of the magnetic tail. The I M F orien- 
tation included in this study is perpendicular to the solar 
wind flow direction. The numerical computations have 
been carried on for a long time which allow us to believe 
that we approached a quasi-steady state. 

The perpendicular case, although similar in some 
aspects, differs in a general behavior to the parallel case. 
In particular, it has been found that the plasma dynamics 
is hardly affected by the strength of the magnetic field 
which originated from the solar wind. The location of the 
bow shock is virtually the same for different values of the 
plasma /?. The IMF lines pile up against the atmospheric 
plasma and drape over the planet to form a magnetoiaii 
with two lobes of opposite magnetic poiaruv separated b\ 
a plasma sheet. The parallel magnetic field causes me 
plasma dynamics to be highly sensitive to a value of the 
plasma /? (Murawski and Steinoifson, 1995). 

The numerical results indicate that the presence of M L 
causes the bow shock to move outward from the planet 
but not so far as observed at Venus during the period 
of maximum solar activity. Some part of the outward 
displacement of the bow shock from that indicated for 
flows without ML thus appears to be accounted for. but 
additional considerations are required to achieve a good 
quantitative agreement with the observational data set 
Whether the main source of these differences arises from 
shortcomings in the model itself or in the interpretation 
of observational data remains uncertain, and is the subject 
of continuing studv. 

It should be recalled in this regard that the mass loading 
model presented in this work concerns only the ions 
pickup due to photoiomzation. The neglect of other ions' 
production processes such as proton-neutral charge ex- 
change and impact ionization could result in unuer- 
estimations of the ion production rate. As a conseuucnce. 
me inclusion of suen additional processes would uncount- 
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